#include "media\shaders\renderconfig.cfg"
#include "FluidGridSimulator.h"
#include <math.h>
#include <IL/il.h>

FluidGridSimulator::FluidGridSimulator(int m, int n, int l, float sizeX, float sizeY, float sizeZ, float posX, float posY, float posZ)
{
	//member initialization
	sprintf(displayModeString,"3D density display");
	lockExternalForces = false;
	renderSliceCount = 64;
	viscosity = 0;
	iterationCount = 50;
	earlyZCull = false;
	addObstacleInfluence = false;
	renderBoundingCube = true;
	type = 0;	
	this ->sizeX = sizeX;
	this ->sizeY = sizeY;
	this ->sizeZ = sizeZ;
	this->posX = posX;
	this->posY = posY;
	this->posZ = posZ;

	m_size[0]=sizeX / 2; m_size[1] = sizeZ / 2; m_size[2] = sizeY / 2;		//size of the bounding box
	m_pos[0]=posX; m_pos[1] = posZ; m_pos[2] = posY;						//position of the bounding box

	externalForces = false;
	pause = false;
	isCamChanged = true;
	computeBuoyancy = true;
	computeDensGrad = false;
	computeVorticity = true;
	realSmoke = true;
	displaymode = 0;
	strokeTexType = 0;
	innerStrokeType = 0;
	nprStyle = 0;
	res_m = m;
	res_n = n;
	res_l = l;
	//no slip condition (value at boundary = 0)
	velocityBoundScale = -1;
	//pure Neumann condition (derivate at boundary = 0)
	pressureBoundScale = 1;
	//these options were only used for testing, 
	//they should always have the following values in the final version:
	clearPressurePerFrame = true;
	computePressureBoundPerIter = false;

	//initial impulse position
	initializeImpulse();
	m_time = 0.f;
	//wang tiles load
	wang.loadTileSet("data/tileset.dat");

	//compute slice resolution
	int M = res_m + 2; //2 pixel bound edges
	int N = res_n + 2; //2 pixel bound edges
	int L = res_l + 2; //2 boundary layers
	if(res_l == 1) L = 1;

	//compute number of grid slice tiles in each row and column of the flat3D textures
	numTilesX = floor(sqrt((double)L));
	numTilesY = ceil((double)(L - numTilesX * numTilesX)/ numTilesX) + numTilesX;

	//Create the inner grid layers
	gridLayers = new GridLayer*[l];
	int rownum = 0;
	int columnum = 1;
	for(int i = 0; i < res_l; i++)
	{
		gridLayers[i] = new GridLayer(res_m, res_n, res_l, i, numTilesX, numTilesY, rownum, columnum);
		columnum++;
		if(columnum == numTilesX)
		{
			columnum = 0;
			rownum ++;
		}
	}
	//create the two boundary layers
	boundLayers[0] = new BoundaryLayer(res_m, res_n, res_l, true, numTilesX, numTilesY);
	boundLayers[1] = new BoundaryLayer(res_m, res_n, res_l, false, numTilesX, numTilesY);

	// final flat3D texture resolution
	int fullResX = M * numTilesX;
	int fullResY = N * numTilesY;

	//initialize flat3D render textures
	velocity = new MultiIteratedRenderTexture(fullResX, fullResY);
	externalVelocity = new RenderTexture(fullResX, fullResY, false);
	density = new MultiIteratedRenderTexture(fullResX, fullResY);
	renderedDensity = new IteratedRenderTexture(fullResX, fullResY, false);
	hessian_r1 = new RenderTexture(fullResX, fullResY, false);
	hessian_r2 = new RenderTexture(fullResX, fullResY, false);
	hessian_r3 = new RenderTexture(fullResX, fullResY, false);
	externalDensity = new RenderTexture(fullResX, fullResY, false);
	pressure = new IteratedRenderTexture(fullResX, fullResY, true);
	velocityDiv = new RenderTexture(fullResX, fullResY, false);
	vorticity = new IteratedRenderTexture(fullResX, fullResY, false);
	vorticityForce = new RenderTexture(fullResX, fullResY, false);
	obstacleTexture = 0;

	// create helper texture
	createLayerOffsetTexture();
	
	//create color band texture
	createColorBand();

	//initialize Cg programs
	s_TexturedProgram = new CgProgram();
	s_TexturedProgram->SetProgramFiles(0,"media//shaders//textured_ps.cg");
	s_TexturedProgram->SetProgramEntries(0,"textured_PS");
	s_TexturedProgram->InitPrograms();
	s_ColoredProgram = new CgProgram(0,0,"media//shaders//ConstantColor.cg", "constantColor_PS");
	s_RenderQuadVertexProgram = new CgProgram("media//shaders//render_quad_vs.cg", "main_VS");	
	s_AdvectProgram = new CgProgram(0,0,"media//shaders//Advect_ps.cg", "advect_PS");
	s_BoundaryLayerProgram = new CgProgram(0,0,"media//shaders//Boundary_ps.cg", "boundaryLayer_PS");
	s_BoundaryEdgeProgram = new CgProgram(0,0,"media//shaders//Boundary_ps.cg", "boundaryEdge_PS");
	s_DivergenceProgram = new CgProgram(0,0,"media//shaders//Divergence_ps.cg", "divergence_PS");
	s_JacobiIterProgram = new CgProgram(0,0,"media//shaders//Jacobi_ps.cg", "Jacobi_PS");
	s_ProjectProgram = new CgProgram(0,0,"media//shaders//Project_ps.cg", "projection_PS");
	s_AddForceProgram = new CgProgram(0,0,"media//shaders//addForce_ps.cg", "addVelocitySpot_PS");
	s_AddDensityProgram = new CgProgram(0,0,"media//shaders//addForce_ps.cg", "addDensitySpot_PS");
	s_VolumeRenderSliceFlat3DProgram = new CgProgram("media//shaders//VolumeRender.cg", "slicesEX_VS","media//shaders//VolumeRender.cg", "slicesFlat3D_PS");
	s_VorticityProgram = new CgProgram(0, 0,"media//shaders//Vorticity_ps.cg", "vorticity_PS");
	s_VorticityForceProgram = new CgProgram(0, 0,"media//shaders//Vorticity_ps.cg", "vorticityForce_PS");
	s_WritePressureAsDepthProgram = new CgProgram(0, 0,"media//shaders//PressureDepth_ps.cg", "writePressureAsDepth");

	//initialize Cg programs
	s_BuoyancyProgram = new CgProgram(0, 0, "media//shaders//addForce_ps.cg", "addBuoyancyForce_PS");
	s_denGraProgram = new CgProgram(0, 0, "media//shaders//DensityGradient_ps.cg", "computeGradient_PS");
	s_denHessianProgram = new CgProgram(0, 0, "media//shaders//DensityGradient_ps.cg", "computeHessian_PS");

	//initialize Cg programs for ray casting
	s_IntersecfrontProgram = new CgProgram("media//shaders//SmokeRender.cg", "intersec_VS", "media//shaders//SmokeRender.cg", "intersec_front_PS");
	s_IntersecbackProgram = new CgProgram("media//shaders//SmokeRender.cg", "intersec_VS", "media//shaders//SmokeRender.cg", "intersec_back_PS");
	s_LocateIntersectProgram = new CgProgram(0,0, "media//shaders//SmokeRender.cg", "locateIntersection_PS");

	//programs for grid point rendering
	s_gridpointProgram = new CgProgram(0, 0, "media//shaders//VolumeRender.cg", "gridpoint_PS");

	//programs for smoothing
	s_smoothProgram = new CgProgram(0, 0, "media//shaders//DensityGradient_ps.cg", "smooth_PS");
	s_smooth2DProgram = new CgProgram(0, 0, "media//shaders//DensityGradient_ps.cg", "smooth2D_PS");

	//programs for rendering npr smoke
	s_realRenderProgram = new CgProgram(0, 0,"media//shaders//VolumeRender.cg", "realRender_PS");
	s_renderTexProgram = new CgProgram(0, 0,"media//shaders//VolumeRender.cg", "renderTex_PS");

	s_alphatestProgram = new CgProgram("media//shaders//VolumeRender.cg", "slices_VS","media//shaders//VolumeRender.cg", "alphaTest_PS");
	s_extractSilhouetteProgram = new CgProgram(0, 0, "media//shaders//VolumeRender.cg", "extractSilhouette_PS");
	s_renderResultTex = new CgProgram(0, 0, "media//shaders//VolumeRender.cg", "renderResultTex_PS");
	//Cg programs for maccormack scheme
	s_advectMacCormack = new CgProgram(0, 0, "media//shaders//Advect_ps.cg", "advectMacCormack_PS");
	//Cg programs for light accumulation
	s_lightProgram = new CgProgram("media//shaders//VolumeRender.cg", "slicesEX_VS","media//shaders//VolumeRender.cg", "lightAccumulate_PS");
	//Cg program for hatching
	s_hatchingProgram = new CgProgram("media//shaders//Hatching.cg", "hatching_VS", "media//shaders//Hatching.cg", "hatching_PS"
		,"media//shaders//Hatching.cg", "hatching_GS");
	s_intertexturesProgram = new CgProgram("media//shaders//PointDis.cg", "intertextures_VS", 
		"media//shaders//PointDis.cg", "intertextures_PS");
	s_edgedetectProgram = new CgProgram(0, 0, "media//shaders//Contour.cg", "edgedetect_PS");
	s_texcontourProgram = new CgProgram("media//shaders/Contour.cg", "texturingcontour_VS", "media//shaders//Contour.cg", "texturingcontour_PS"
		,"media//shaders//Contour.cg", "texturingcontour_GS");
	s_texcontourVProgram = new CgProgram("media//shaders/Contour.cg", "texturingcontour_VS", "media//shaders//Contour.cg", "texturingcontour_PS"
		,"media//shaders//Contour.cg", "texturingcontourV_GS");
	s_toonshadeProgram = new CgProgram(0, 0, "media//shaders//Contour.cg", "toonshade_PS");
	s_NPRcolorizeProgram = new CgProgram("media//shaders//Contour.cg", "NPRcolorize_VS", 
		"media//shaders//Contour.cg", "NPRcolorize_PS", "media//shaders//Contour.cg", "NPRcolorize_GS");
	s_advectSampPoint = new CgProgram("media//shaders//PointDis.cg", "advectSampPoint_VS", 
		"media//shaders//PointDis.cg", "advectSampPoint_PS", 
		"media//shaders//PointDis.cg", "advectSampPoint_GS");
	s_distPointProgram = new CgProgram("media//shaders//PointDis.cg", "distributePoint_VS", 
		"media//shaders//PointDis.cg", "distributePoint_PS");
	s_distPointProgram2v = new CgProgram("media//shaders//PointDis.cg", "distributePoint2v_VS", 
		"media//shaders//PointDis.cg", "distributePoint2v_PS");
	s_markRegionProgram = new CgProgram("media//shaders//PointDis.cg", "markRegion_VS", 
		"media//shaders//PointDis.cg", "markRegion_PS"
		,"media//shaders//PointDis.cg", "markRegion_GS");
	s_stippleProgram = new CgProgram("media//shaders//Stippling.cg", "stippling_VS", 
		"media//shaders//Stippling.cg", "stippling_PS");

	s_rotateSampPoint = new CgProgram("media//shaders//PointDis.cg", "rotateSampPoint_VS", 
		"media//shaders//PointDis.cg", "advectSampPoint_PS");

	s_billboardProgram = new CgProgram("media//shaders//BillBoard.cg", "billboard_VS",
		"media//shaders//BillBoard.cg", "billboard_PS");
	s_solidcontourProgram = new CgProgram("media//shaders//BillBoard.cg", "solidcontour_VS", 
		"media//shaders//BillBoard.cg", "solidcontour_PS"
		,"media//shaders//BillBoard.cg", "solidcontour_GS");

	//initialize flat3D textures
	intersection = NULL;
	buoyancyForce = new RenderTexture(fullResX, fullResY, false);
	temperature = new IteratedRenderTexture(fullResX, fullResY, false);
	externalTemp = new RenderTexture(fullResX, fullResY, false);
	phi_n_hat = new RenderTexture(fullResX, fullResY, false);
	phi_n_1_hat = new RenderTexture(fullResX, fullResY, false);
	lw = lh = 100.f;
	lightbuffer = new RenderTexture(lw, lh, false);

	//create billboard texture
	createPuffBillBoard();

	footprint = NULL;
	resulttex = NULL;
	imgbuff = NULL;
	marktex = NULL;
	ndotetex = NULL;
	edgetex = NULL;
	depthtex = NULL;
	pointCollection = NULL;
	initpointtex = NULL;
	advectCollection = NULL;
	tmpadvecttex = NULL;
	advectCollection_2v = NULL;
	pointCollection_2v = NULL;
	tmpadvecttex_2v = NULL;

	m_disthres = 0.05;//0.01;
	dsampleX = 20;
	dsampleY = 20;
	m_contourtexsize = 20;
	m_sampleRatio = 10;
	nSmooth = 1;
	recsize = 5;
	pointdisalg = 0;
	edgemode = 0;
	joinmode = 0;
	gausssig = 1;
	hatchmode = 0;
	detectMode = 0;
	densityThres = 0.003f;
	intensityThres = 0.05f;
	m_contourthres = 0.01f;
	m_sclstep = 0;

	//get the current ModelViewProj matrix
	glGetDoublev(GL_MODELVIEW_MATRIX, preModelViewProj);
	glMatrixMode(GL_PROJECTION); glPushMatrix();
	glMultMatrixd(preModelViewProj);
	glGetDoublev(GL_PROJECTION_MATRIX, preModelViewProj);
	glMatrixMode(GL_PROJECTION); glPopMatrix();

	dt = 0.003;//0.003;//0.007;
	//generate k first density field
	int textureNumber = reinterpret_cast<MultiIteratedRenderTexture*>(density)->getTextureNumber();
	for (int i = 0; i < textureNumber - 1; i++)
	{
		setVelocityBoundScale(-1);
		setPressureBoundScale(1);
		update();
	}

}

FluidGridSimulator::~FluidGridSimulator(void)
{
	glDeleteTextures(1, &layerOffsetTexture);
	glDeleteTextures(1, &strokeTexture);
	glDeleteTextures(1, &nprcolorBand);
	glDeleteTextures(1, &halftoneBand);
	glDeleteTextures(1, &splatTexture);
	glDeleteTextures(1, &opaquediskTexture);
	glDeleteTextures(1, &contourstrokeTexture);

	delete boundLayers[0];
	delete boundLayers[1];
	for(int i = 0; i < res_l; i++)
		delete gridLayers[i];
	delete[] gridLayers;

	delete velocity;
	delete externalVelocity;
	delete density;
	delete externalDensity;
	delete pressure;
	delete velocityDiv;
	delete vorticity;
	delete vorticityForce;

	delete s_TexturedProgram;
	delete s_ColoredProgram;
	delete s_RenderQuadVertexProgram;	
	delete s_AdvectProgram;
	delete s_BoundaryLayerProgram;
	delete s_BoundaryEdgeProgram;
	delete s_DivergenceProgram;
	delete s_JacobiIterProgram;
	delete s_ProjectProgram;
	delete s_AddForceProgram;
	delete s_AddDensityProgram;
	delete s_VolumeRenderSliceFlat3DProgram;
	delete s_VorticityProgram;
	delete s_VorticityForceProgram;
	delete s_WritePressureAsDepthProgram;

	delete s_BuoyancyProgram;
	delete s_denGraProgram;
	delete s_denHessianProgram;

	delete s_IntersecfrontProgram;
	delete s_IntersecbackProgram;
	delete s_LocateIntersectProgram;
	delete s_gridpointProgram;
	delete s_smoothProgram;
	delete s_smooth2DProgram;
	delete s_realRenderProgram;
	delete s_renderTexProgram;
	delete s_alphatestProgram;
	delete s_extractSilhouetteProgram;
	delete s_renderResultTex;
	delete s_intertexturesProgram;
	delete s_edgedetectProgram;
	delete s_texcontourProgram;
	delete s_texcontourVProgram;
	delete s_NPRcolorizeProgram;
	delete s_toonshadeProgram;
	delete s_advectSampPoint;
	delete s_distPointProgram;
	delete s_markRegionProgram;
	delete s_stippleProgram;
	delete s_hatchingProgram;
	delete s_rotateSampPoint;
	delete s_billboardProgram;
	delete s_solidcontourProgram;
	delete s_distPointProgram2v;
	//delete s_advectMacCormack;

	delete buoyancyForce;
	delete externalDensity;
	delete externalTemp;
	delete temperature;
	//delete phi_n_1_hat;
	//delete phi_n_hat;
	delete footprint;
	delete resulttex;
	delete imgbuff;
	delete marktex;
	delete ndotetex;
	delete edgetex;
	delete depthtex;
	delete pointCollection;
	delete initpointtex;
	delete billBoard;
	delete advectCollection;
	delete tmpadvecttex;
	delete advectCollection_2v;
	delete pointCollection_2v;
	delete tmpadvecttex_2v;

	glDeleteLists(dispListId, 1);
	glDeleteLists(samplingListId, 1);
	
	delete[] randIntensity;
}

// set external sources to zero
void FluidGridSimulator::clearExternalForces()
{
	externalVelocity->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);
	externalVelocity->endRender();
	externalDensity->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);
	externalDensity->endRender();

	externalTemp->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);
	externalTemp->endRender();
}

// simulation of timestep dt
void FluidGridSimulator::update()
{	
	s_RenderQuadVertexProgram->Enable();


	//advect density and velocity	
	Advect(dt, velocity, velocity, velocityBoundScale);
	Advect(dt, density, velocity, 0);
	//Advect(dt, temperature, velocity, 0);

	// add external forces if needed
	addGaussianBall();
	if(externalForces)
		AddExternalForces();

	//add buoyancy force
	if(computeBuoyancy)
		addBuoyancyForces(dt);

	// add obstacle influence if needed
	if(addObstacleInfluence)
		AddObstacleInfluence();

	// compute and add vorticity influence if needed
	if(computeVorticity)
		ComputeVorticity(dt);	

	// diffuse velocity if needed
	if(viscosity)
		DiffuseVelocity(dt);

	// compute pressure and do the final projection step
	ComputePressure();
	Project();

	s_RenderQuadVertexProgram->Disable();
}

void FluidGridSimulator::ComputeBoundary(IteratedRenderTexture* media, float borderScale, bool updateUpperBound, bool updateLowerBound)
{	
	//set shader programs and parameters for boundary layer rendering
	s_BoundaryLayerProgram->SetFragmentTexParam("media", media->getTextureID());
	s_BoundaryLayerProgram->Enable();
	s_BoundaryLayerProgram->EnableFragmentTexParam("media");
	s_BoundaryLayerProgram->SetFragmentParam1f("scale", borderScale);
	//render 2 boundary layers
	if(updateUpperBound)
		boundLayers[0]->DrawQuad();
	if(updateLowerBound)
		boundLayers[1]->DrawQuad();
	s_BoundaryLayerProgram->DisableFragmentTexParam("media");
	s_BoundaryLayerProgram->Disable();
	//set shader programs and parameters for boundary edge rendering
	s_BoundaryEdgeProgram->SetFragmentTexParam("media", media->getTextureID());
	s_BoundaryEdgeProgram->Enable();
	s_BoundaryEdgeProgram->EnableFragmentTexParam("media");
	s_BoundaryEdgeProgram->SetFragmentParam1f("scale", borderScale);
	for(int i = 0; i < res_l; i++)
		gridLayers[i]->DrawBounds();
	s_BoundaryEdgeProgram->DisableFragmentTexParam("media");
	s_BoundaryEdgeProgram->Disable();		
}

//used to advect velocity and density
void FluidGridSimulator::Advect(float dt, IteratedRenderTexture* media, IteratedRenderTexture* velocity, float borderScale)
{
	bool updateUpper = true;
	bool updateLower = true;

	//distributed simulation
	if(type == 1)
		updateUpper = false;
	if(type == 2)
	{
		updateUpper = false;
		updateLower = false;
	}
	if(type == 3)
		updateLower = false;	

	//media->startRender();
	////compute boundary values
	//ComputeBoundary(media, borderScale, updateUpper, updateLower);

	////set shader programs and parameters for inner layer rendering
	//s_AdvectProgram->SetFragmentTexParam("velocity", velocity->getTextureID());
	//s_AdvectProgram->SetFragmentTexParam("media", media->getTextureID());
	//s_AdvectProgram->SetFragmentTexParam("offsetTexture", layerOffsetTexture);
	//s_AdvectProgram->Enable();
	//s_AdvectProgram->EnableFragmentTexParam("velocity");
	//s_AdvectProgram->EnableFragmentTexParam("media");
	//s_AdvectProgram->EnableFragmentTexParam("offsetTexture");
	//s_AdvectProgram->SetFragmentParam1f("dt", dt);
	//s_AdvectProgram->SetFragmentParam2f("layerResolution", res_m, res_n);

	////render inner grid layers
	//for(int i = 0; i < res_l; i++)
	//{
	//	s_AdvectProgram->SetFragmentParam1f("l", i + 1);
	//	gridLayers[i]->DrawQuad();
	//}

	//s_AdvectProgram->DisableFragmentTexParam("velocity");
	//s_AdvectProgram->DisableFragmentTexParam("media");
	//s_AdvectProgram->DisableFragmentTexParam("offsetTexture");
	//s_AdvectProgram->Disable();

	//media->endRender();

	//////////////////////////////////////////
	//MacCormack advection scheme
	
	//compute phi n+1 hat
	phi_n_1_hat->startRender();
	//renderTexturedQuad(media);
	s_AdvectProgram->SetFragmentTexParam("velocity", velocity->getTextureID());
	s_AdvectProgram->SetFragmentTexParam("media", media->getTextureID());
	s_AdvectProgram->SetFragmentTexParam("offsetTexture", layerOffsetTexture);
	s_AdvectProgram->Enable();
	s_AdvectProgram->EnableFragmentTexParam("velocity");
	s_AdvectProgram->EnableFragmentTexParam("media");
	s_AdvectProgram->EnableFragmentTexParam("offsetTexture");
	s_AdvectProgram->SetFragmentParam1f("dt", dt);
	s_AdvectProgram->SetFragmentParam2f("layerResolution", res_m, res_n);
	for(int i = 0; i < res_l; i++)
	{
		s_AdvectProgram->SetFragmentParam1f("l", i + 1);
		gridLayers[i]->DrawQuad();
	}
	s_AdvectProgram->DisableFragmentTexParam("velocity");
	s_AdvectProgram->DisableFragmentTexParam("media");
	s_AdvectProgram->DisableFragmentTexParam("offsetTexture");
	s_AdvectProgram->Disable();
	phi_n_1_hat->endRender();

	//compute phi n hat
	phi_n_hat->startRender();
	s_AdvectProgram->SetFragmentTexParam("velocity", velocity->getTextureID());
	s_AdvectProgram->SetFragmentTexParam("media", phi_n_1_hat->getTextureID());
	s_AdvectProgram->SetFragmentTexParam("offsetTexture", layerOffsetTexture);
	s_AdvectProgram->Enable();
	s_AdvectProgram->EnableFragmentTexParam("velocity");
	s_AdvectProgram->EnableFragmentTexParam("media");
	s_AdvectProgram->EnableFragmentTexParam("offsetTexture");
	s_AdvectProgram->SetFragmentParam1f("dt", -dt);
	s_AdvectProgram->SetFragmentParam2f("layerResolution", res_m, res_n);
	for(int i = 0; i < res_l; i++)
	{
		s_AdvectProgram->SetFragmentParam1f("l", i + 1);
		gridLayers[i]->DrawQuad();
	}
	s_AdvectProgram->DisableFragmentTexParam("velocity");
	s_AdvectProgram->DisableFragmentTexParam("media");
	s_AdvectProgram->DisableFragmentTexParam("offsetTexture");
	s_AdvectProgram->Disable();
	phi_n_hat->endRender();

	//complete advection
	media->startRender();
	s_advectMacCormack->SetFragmentTexParam("velocity", velocity->getTextureID());
	s_advectMacCormack->SetFragmentTexParam("phi_n", media->getTextureID());
	s_advectMacCormack->SetFragmentTexParam("phi_n_1_hat", phi_n_1_hat->getTextureID());
	s_advectMacCormack->SetFragmentTexParam("phi_n_hat", phi_n_hat->getTextureID());
	s_advectMacCormack->SetFragmentTexParam("offsetTexture", layerOffsetTexture);
	s_advectMacCormack->Enable();
	s_advectMacCormack->EnableFragmentTexParam("velocity");
	s_advectMacCormack->EnableFragmentTexParam("phi_n");
	s_advectMacCormack->EnableFragmentTexParam("phi_n_1_hat");
	s_advectMacCormack->EnableFragmentTexParam("phi_n_hat");
	s_advectMacCormack->EnableFragmentTexParam("offsetTexture");
	s_advectMacCormack->SetFragmentParam1f("dt", dt);
	s_advectMacCormack->SetFragmentParam2f("layerResolution", res_m, res_n);
	for(int i = 0; i < res_l; i++)
	{
		s_advectMacCormack->SetFragmentParam1f("l", i + 1);
		gridLayers[i]->DrawQuad();
	}
	s_advectMacCormack->DisableFragmentTexParam("velocity");
	s_advectMacCormack->DisableFragmentTexParam("phi_n");
	s_advectMacCormack->DisableFragmentTexParam("phi_n_1_hat");
	s_advectMacCormack->DisableFragmentTexParam("phi_n_hat");
	s_advectMacCormack->DisableFragmentTexParam("offsetTexture");
	media->endRender();
}

// Add external sources to the velocity and density fields
// Implemented with fullscreen texture quad rendering and additive blending
void FluidGridSimulator::AddExternalForces()
{	
	glEnable(GL_BLEND);
	glBlendFunc(GL_ONE, GL_ONE);

	velocity->startCopy();
	renderTexturedQuad(externalVelocity);
	velocity->endCopy();

	density->startCopy();
	renderTexturedQuad(externalDensity);
	density->endCopy();

	temperature->startCopy();
	renderTexturedQuad(externalTemp);
	temperature->endCopy();

	glDisable(GL_BLEND);

}

// Add obstacle influance to the density and velocity fields
// implemented with fullscreen quad rendering and alpha blending
void FluidGridSimulator::AddObstacleInfluence()
{
	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
	//glBlendFunc(GL_ONE, GL_ONE);	

	velocity->startCopy();
	renderTexturedQuad(obstacleTexture);
	velocity->endCopy();

	density->startCopy();
	renderTexturedQuad(obstacleTexture);
	density->endCopy();

	glDisable(GL_BLEND);	
}

// add a uniform (equal in the whole grid) density and/or density to the external forces
void FluidGridSimulator::addExternalUniformSource(float density, float temp, float velocityX, float velocityY, float velocityZ)
{
	externalForces = true;

	s_RenderQuadVertexProgram->Enable();

	if(density > 0)
	{
		externalDensity->startRender();

		glEnable(GL_BLEND);
		glBlendFunc(GL_ONE, GL_ONE);

		s_ColoredProgram->Enable();
		s_ColoredProgram->SetFragmentParam4f("color", density, density, density, 0);

		this->renderFullscreenQuad();

		s_ColoredProgram->Disable();

		glDisable(GL_BLEND);

		externalDensity->endRender();
	}
	if( (velocityX * velocityX + velocityY * velocityY + velocityZ * velocityZ) > 0 )
	{
		externalVelocity->startRender();

		glEnable(GL_BLEND);
		glBlendFunc(GL_ONE, GL_ONE);

		s_ColoredProgram->Enable();
		s_ColoredProgram->SetFragmentParam4f("color", velocityX, velocityY, velocityZ, 0);

		this->renderFullscreenQuad();

		s_ColoredProgram->Disable();

		glDisable(GL_BLEND);

		externalVelocity->endRender();
	}

	if (temp > 0)
	{
		externalTemp->startRender();

		glEnable(GL_BLEND);
		glBlendFunc(GL_ONE, GL_ONE);

		s_ColoredProgram->Enable();
		s_ColoredProgram->SetFragmentParam4f("color", temp, temp, temp, 0);

		this->renderFullscreenQuad();

		s_ColoredProgram->Disable();

		glDisable(GL_BLEND);

		externalTemp->endRender();
	}

	s_RenderQuadVertexProgram->Disable();
}
// add buoyancy force
void FluidGridSimulator::addBuoyancyForces(float dt)
{
	bool updateUpper = true;
	bool updateLower = true;
	if(type == 1)
		updateUpper = false;
	if(type == 2)
	{
		updateUpper = false;
		updateLower = false;
	}
	if(type == 3)
		updateLower = false;

	//copy inner pixels and compute boundary pixels
	velocity->startRender();
	renderTexturedQuad(velocity);
	ComputeBoundary(velocity, velocityBoundScale, updateUpper, updateLower);	
	velocity->endRender();

	//random direction
	Vector3f forcedir;
	forcedir.x = (float)rand()/RAND_MAX - 0.5f;
	forcedir.y = (float)rand()/RAND_MAX - 0.5f;
	forcedir.z = (float)rand()/RAND_MAX - 0.5f;
	forcedir.normalize();

	// compute force caused by buoyancy
	buoyancyForce->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);
	s_BuoyancyProgram->SetFragmentTexParam("density", density->getTextureID());
	s_BuoyancyProgram->SetFragmentTexParam("temperature", temperature->getTextureID());
	s_BuoyancyProgram->SetFragmentParam3f("forcedir", forcedir.x, forcedir.y, forcedir.z);
	s_BuoyancyProgram->Enable();
	s_BuoyancyProgram->EnableFragmentTexParam("density");
	s_BuoyancyProgram->EnableFragmentTexParam("temperature");
	s_BuoyancyProgram->SetFragmentParam1f("dt", dt);
	for(int i = 0; i < res_l; i++)
		gridLayers[i]->DrawQuad();
	s_BuoyancyProgram->DisableFragmentTexParam("density");
	s_BuoyancyProgram->DisableFragmentTexParam("temperature");
	s_BuoyancyProgram->Disable();
	buoyancyForce->endRender();

	// add the vorticity force to the velocity field
	velocity->startCopy();
	glEnable(GL_BLEND);
	glBlendFunc(GL_ONE, GL_ONE);
	renderTexturedQuad(buoyancyForce);
	glDisable(GL_BLEND);
	velocity->endCopy();
}

// add a point source to the external sources
// implemented with fullscreen quad rendering and additive blending
static float lilrand()
{
	return (rand()/float(RAND_MAX) - 0.5f)*5.0f;
}

void FluidGridSimulator::addGaussianBall()
{
	externalForces = true;
	float impulseSize = 2;
	float saturation = 0.78f;

	Vector3f dir;
	dir.x = 0.f; dir.y = 0.f; dir.z = -1.f;
	dir.normalize();

	//m_impulse += dir * 1.f;
	//m_impulse.y = res_m / 6.f;

	//add density
	float density = 0.3f;//0.05f;
	externalDensity->startRender();
	glClearColor(0, 0, 0, 0);
	glClear(GL_COLOR_BUFFER_BIT);
	s_AddDensityProgram->Enable();
	s_AddDensityProgram->SetFragmentParam3f("center", m_impulse.x + lilrand(), 
		m_impulse.y + lilrand(), m_impulse.z + lilrand());
	s_AddDensityProgram->SetFragmentParam1f("size", impulseSize);
	s_AddDensityProgram->SetFragmentParam4f("splatColor", density, density, density, 0.f);
	for(int i = 0; i < res_l; i++)
	{
		s_AddDensityProgram->SetFragmentParam1f("l", i + 1);
		gridLayers[i]->DrawQuad();
	}
	s_AddDensityProgram->Disable();
	externalDensity->endRender();

	//add velocity
	impulseSize = 2.f;
	externalVelocity->startRender();
	glClearColor(0, 0, 0, 0);
	glClear(GL_COLOR_BUFFER_BIT);
	//render points of forces
	s_AddDensityProgram->Enable();
	s_AddDensityProgram->SetFragmentParam3f("center", m_impulse.x + lilrand(), m_impulse.y + lilrand(), m_impulse.z + lilrand());
	s_AddDensityProgram->SetFragmentParam1f("size", impulseSize);
	float force = 500.f;// + 200.f*(float)rand()/RAND_MAX;//500.f;
	s_AddDensityProgram->SetFragmentParam4f("splatColor", -force*dir.x, -force*dir.y, -force*dir.z, 0.f);
	for(int i = 0; i < res_l; i++)
	{
		s_AddDensityProgram->SetFragmentParam1f("l", i + 1);
		gridLayers[i]->DrawQuad();
	}
	s_AddDensityProgram->Disable();
	externalVelocity->endRender();

	//add controllable wind
	if (m_isWind)
	{
		impulseSize = 5;
		externalVelocity->startRender();
		glEnable(GL_BLEND);	glBlendFunc(GL_ONE, GL_ONE);
		s_AddDensityProgram->Enable();
		s_AddDensityProgram->SetFragmentParam3f("center", m_windsource[0], m_windsource[1], m_windsource[2]);
		s_AddDensityProgram->SetFragmentParam1f("size", impulseSize);
		force = 50.f;
		s_AddDensityProgram->SetFragmentParam4f("splatColor", force*m_winddir[0], force*m_winddir[1], force*m_winddir[2] , 0.f);
		for(int i = 0; i < res_l; i++)
		{
			s_AddDensityProgram->SetFragmentParam1f("l", i + 1);
			gridLayers[i]->DrawQuad();
		}
		s_AddDensityProgram->Disable();
		externalVelocity->endRender();
		glDisable(GL_BLEND);
		m_isWind = false;
	}
}

void FluidGridSimulator::addExternalSpotSource(float* position,
											   float* direction,
											   float velocityRadius,
											   float densityRadius,
											   float tempRadius,
											   float velocityScale,
											   float densityScale,
											   float tempScale)
{
	externalForces = true;

	s_RenderQuadVertexProgram->Enable();

	if(velocityScale * velocityRadius)
	{

		//render forces to velocity texture
		//set shader programs and parameters
		externalVelocity->startRender();
		glEnable(GL_BLEND);
		glBlendFunc(GL_ONE, GL_ONE);
		//render points of forces
		s_AddForceProgram->Enable();
		s_AddForceProgram->SetFragmentParam3f("sourcePos", position[0], position[1], position[2]);
		s_AddForceProgram->SetFragmentParam3f("sourceDir", direction[0], direction[1], direction[2]);
		s_AddForceProgram->SetFragmentParam1f("sourceRadius", velocityRadius);
		s_AddForceProgram->SetFragmentParam1f("scale", velocityScale);
		for(int i = 0; i < res_l; i++)
		{
			s_AddForceProgram->SetFragmentParam1f("l", i + 1);
			gridLayers[i]->DrawQuad();
		}
		glDisable(GL_BLEND);
		s_AddForceProgram->Disable();
		externalVelocity->endRender();
	}

	if(densityScale * tempRadius)
	{	
		externalDensity->startRender();
		glEnable(GL_BLEND);
		glBlendFunc(GL_ONE, GL_ONE);

		s_AddDensityProgram->Enable();
		s_AddDensityProgram->SetFragmentParam3f("sourcePos", position[0], position[1], position[2]);
		s_AddDensityProgram->SetFragmentParam1f("sourceRadius", tempRadius);
		s_AddDensityProgram->SetFragmentParam1f("scale", densityScale);
		for(int i = 0; i < res_l; i++)
		{
			s_AddDensityProgram->SetFragmentParam1f("l", i + 1);
			gridLayers[i]->DrawQuad();
		}
		glDisable(GL_BLEND);
		s_AddDensityProgram->Disable();
		externalDensity->endRender();
	}

	if(tempScale * densityRadius)
	{	
		externalTemp->startRender();
		glEnable(GL_BLEND);
		glBlendFunc(GL_ONE, GL_ONE);

		s_AddDensityProgram->Enable();
		s_AddDensityProgram->SetFragmentParam3f("sourcePos", position[0], position[1], position[2]);
		s_AddDensityProgram->SetFragmentParam1f("sourceRadius", densityRadius);
		s_AddDensityProgram->SetFragmentParam1f("scale", tempScale);
		for(int i = 0; i < res_l; i++)
		{
			s_AddDensityProgram->SetFragmentParam1f("l", i + 1);
			gridLayers[i]->DrawQuad();
		}
		glDisable(GL_BLEND);
		s_AddDensityProgram->Disable();
		externalTemp->endRender();
	}

	s_RenderQuadVertexProgram->Disable();
}

void FluidGridSimulator::ComputeVorticity(float dt)
{
	bool updateUpper = true;
	bool updateLower = true;
	if(type == 1)
		updateUpper = false;
	if(type == 2)
	{
		updateUpper = false;
		updateLower = false;
	}
	if(type == 3)
		updateLower = false;

	//copy inner pixels and compute boundary pixels
	velocity->startRender();
	renderTexturedQuad(velocity);
	ComputeBoundary(velocity, velocityBoundScale, updateUpper, updateLower);	
	velocity->endRender();

	// compute the vorticity of the velocity field
	vorticity->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);
	s_VorticityProgram->SetFragmentTexParam("velocity", velocity->getTextureID());
	s_VorticityProgram->Enable();
	s_VorticityProgram->EnableFragmentTexParam("velocity");
	for(int i = 0; i < res_l; i++)
		gridLayers[i]->DrawQuad();				
	s_VorticityProgram->DisableFragmentTexParam("velocity");
	s_VorticityProgram->Disable();
	vorticity->endRender();

	// compute velocity caused by the vorticity field
	vorticityForce->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);
	s_VorticityForceProgram->SetFragmentTexParam("vorticity", vorticity->getTextureID());
	s_VorticityForceProgram->Enable();
	s_VorticityForceProgram->EnableFragmentTexParam("vorticity");
	s_VorticityForceProgram->SetFragmentParam1f("dt", dt);
	for(int i = 0; i < res_l; i++)
		gridLayers[i]->DrawQuad();
	s_VorticityForceProgram->DisableFragmentTexParam("vorticity");
	s_VorticityForceProgram->Disable();
	vorticityForce->endRender();

	// add the vorticity force to the velocity field
	velocity->startCopy();
	glEnable(GL_BLEND);
	glBlendFunc(GL_ONE, GL_ONE);
	renderTexturedQuad(vorticityForce);
	glDisable(GL_BLEND);
	velocity->endCopy();
}

void FluidGridSimulator::DiffuseVelocity(float dt)
{
	float alpha = 1.0 / (viscosity * dt * 1.0);
	//Jacoby iteration
	for(int i = 0; i< iterationCount; i++)
	{
		velocity->startRender();

		if(computePressureBoundPerIter)
			ComputeBoundary(velocity, velocityBoundScale);

		s_JacobiIterProgram->SetFragmentTexParam("x", velocity->getTextureID());
		s_JacobiIterProgram->SetFragmentTexParam("b", velocity->getTextureID());
		s_JacobiIterProgram->Enable();
		s_JacobiIterProgram->EnableFragmentTexParam("x");
		s_JacobiIterProgram->EnableFragmentTexParam("b");
		s_JacobiIterProgram->SetFragmentParam1f("alpha", alpha);
		s_JacobiIterProgram->SetFragmentParam1f("beta", 6 + alpha);
		for(int l = 0; l < res_l; l++)
			gridLayers[l]->DrawQuad();
		s_JacobiIterProgram->DisableFragmentTexParam("x");
		s_JacobiIterProgram->DisableFragmentTexParam("b");
		s_JacobiIterProgram->Disable();		
		velocity->endRender();
	}

}

void FluidGridSimulator::ComputePressure()
{	
	bool updateUpper = true;
	bool updateLower = true;
	if(type == 1)
		updateUpper = false;
	if(type == 2)
	{
		updateUpper = false;
		updateLower = false;
	}
	if(type == 3)
		updateLower = false;

	//compute velocity divergence
	velocityDiv->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);
	s_DivergenceProgram->SetFragmentTexParam("velocity", velocity->getTextureID());
	s_DivergenceProgram->Enable();
	s_DivergenceProgram->EnableFragmentTexParam("velocity");

	for(int i = 0; i < res_l; i++)
		gridLayers[i]->DrawQuad();

	s_DivergenceProgram->DisableFragmentTexParam("velocity");
	s_DivergenceProgram->Disable();
	velocityDiv->endRender();

	// Early Z culling if needed
	if(earlyZCull)
	{
		//write pressure values to the depth buffer
		pressure->startRender();	
		glClearDepth(1.0);
		glClear(GL_DEPTH_BUFFER_BIT);
		glEnable(GL_DEPTH_TEST);
		renderTexturedQuad( pressure, 1, 0, s_WritePressureAsDepthProgram);
		pressure->endRender();
	}
	glDisable(GL_DEPTH_TEST);

	//zero the pressure
	if(clearPressurePerFrame)
	{
		pressure->startRender();
		//clear whole texture to zero
		glClearColor(0,0,0,0);
		glClear(GL_COLOR_BUFFER_BIT);
		//mark obstacle object voxels (alpha = 1: inside object, alpha = 0: outside object)
		if(addObstacleInfluence)
		{
			glColorMask(false,false,false,true);
			renderTexturedQuad(velocity);
			glColorMask(true,true,true,true);
		}
		pressure->endRender();
	}


	//disable depth write
	if(earlyZCull)
		glEnable(GL_DEPTH_TEST);
	glDepthMask(false);
	//do Jacobi iterations
	for(int i = 0; i < iterationCount; i++)
	{		
		pressure->startRender();
		//only for testing, always false
		if(computePressureBoundPerIter)
			ComputeBoundary(pressure, pressureBoundScale, updateUpper, updateLower);

		s_JacobiIterProgram->SetFragmentTexParam("x", pressure->getTextureID());
		s_JacobiIterProgram->SetFragmentTexParam("b", velocityDiv->getTextureID());
		s_JacobiIterProgram->Enable();
		s_JacobiIterProgram->EnableFragmentTexParam("x");
		s_JacobiIterProgram->EnableFragmentTexParam("b");
		s_JacobiIterProgram->SetFragmentParam1f("alpha", -1);
		s_JacobiIterProgram->SetFragmentParam1f("beta", 6);

		//depth value is increased with the iteration count to cull out pixels with higher and higher pressure values
		float d = (float)i / (float)iterationCount * 2.0 - 1.0;			
		for(int l = 0; l < res_l; l++)
			gridLayers[l]->DrawQuad(d);
		s_JacobiIterProgram->DisableFragmentTexParam("x");
		s_JacobiIterProgram->DisableFragmentTexParam("b");
		s_JacobiIterProgram->Disable();	

		pressure->endRender();
	}
	glDepthMask(true);
	glDisable(GL_DEPTH_TEST);

	//always execute, compute boundary
	if(!computePressureBoundPerIter)
	{
		pressure->startRender();
		//			renderTexturedQuad( pressure, 1, 0);
		ComputeBoundary(pressure, pressureBoundScale, updateUpper, updateLower);
		pressure->endRender();
	}
}

// final projection: substract gradient of the pressure field from the velocity field to get a divergence free velocity field
void FluidGridSimulator::Project()
{
	velocity->startRender();
	s_ProjectProgram->SetFragmentTexParam("w", velocity->getTextureID());
	s_ProjectProgram->SetFragmentTexParam("p", pressure->getTextureID());
	s_ProjectProgram->Enable();
	s_ProjectProgram->EnableFragmentTexParam("w");
	s_ProjectProgram->EnableFragmentTexParam("p");
	for(int l = 0; l < res_l; l++)
		gridLayers[l]->DrawQuad();
	s_ProjectProgram->DisableFragmentTexParam("w");
	s_ProjectProgram->DisableFragmentTexParam("p");	
	s_ProjectProgram->Disable();
	velocity->endRender();

	velocity->startRender();
	renderTexturedQuad(velocity);
	ComputeBoundary(velocity, velocityBoundScale);
	velocity->endRender();
}

// renders a fullscreen quad
void FluidGridSimulator::renderFullscreenQuad()
{
	float dimX = (res_m + 2) * numTilesX;
	float dimY = (res_n + 2) * numTilesY;

	glBegin(GL_QUADS);
	glTexCoord2f(0,0);
	glVertex3f(-1,-1,0);
	glTexCoord2f(dimX,0);
	glVertex3f(1,-1,0);
	glTexCoord2f(dimX,dimY);
	glVertex3f(1,1,0);
	glTexCoord2f(0,dimY);
	glVertex3f(-1,1,0);		
	glEnd();
}

// renders a fullscreen textured quad
// colorTex : texture
// trScale  : texture value scale
// trTrans  : texture value offset
// program  : Cg program to use; if 0 default texturing program is used
void FluidGridSimulator::renderTexturedQuad(IRenderTarget* colorTex, float trScale, float trTrans, CgProgram* program)
{	
	if(program == 0)
		program = s_TexturedProgram;

	program->SetFragmentTexParam("colorMap", colorTex->getTextureID());
	program->Enable();
	program->EnableFragmentTexParam("colorMap");
	program->SetFragmentParam2f("transform", trScale, trTrans);
	float dimX = (res_m + 2) * numTilesX;
	float dimY = (res_n + 2) * numTilesY;

	renderFullscreenQuad();

	program->DisableFragmentTexParam("colorMap");
	program->Disable();	
}

//smoothing
void FluidGridSimulator::smoothingSmokeData()
{
	//smoothing
	GLuint* texlist = reinterpret_cast<MultiIteratedRenderTexture*>(density)->getTextureIdList();
	GLuint* vellist = reinterpret_cast<MultiIteratedRenderTexture*>(velocity)->getTextureIdList();
	s_RenderQuadVertexProgram->Enable();
	renderedDensity->startRender();
	s_smoothProgram->SetFragmentTexParam("image_l1", texlist[0]);	//input density data as 2D texture
	s_smoothProgram->SetFragmentTexParam("image_c", texlist[2]);
	s_smoothProgram->SetFragmentTexParam("image_r1", texlist[1]);
	//s_smoothProgram->SetFragmentTexParam("velocity_c", vellist[2]);
	//s_smoothProgram->SetFragmentTexParam("velocity_l1", vellist[0]);
	s_smoothProgram->SetFragmentTexParam("offsetTexture", layerOffsetTexture);
	s_smoothProgram->SetFragmentParam1f("type", 1);		//dilation
	s_smoothProgram->Enable();
	for (int l = 0; l < res_l; l++)
	{
		s_smoothProgram->SetFragmentParam1f("l", l+1);
		gridLayers[l]->DrawQuad();
	}
	renderedDensity->endRender();
	s_smoothProgram->Disable();
	s_RenderQuadVertexProgram->Disable();

	//smooth more
	for (int i = 0; i < nSmooth; i++)
	{
		s_RenderQuadVertexProgram->Enable();
		renderedDensity->startRender();
		s_smoothProgram->SetFragmentTexParam("image_l1", texlist[0]);	//input density data as 2D texture
		s_smoothProgram->SetFragmentTexParam("image_c", renderedDensity->getTextureID());
		s_smoothProgram->SetFragmentTexParam("image_r1", texlist[1]);
		//s_smoothProgram->SetFragmentTexParam("velocity_c", vellist[2]);
		//s_smoothProgram->SetFragmentTexParam("velocity_l1", vellist[0]);
		s_smoothProgram->SetFragmentTexParam("offsetTexture", layerOffsetTexture);
		s_smoothProgram->SetFragmentParam1f("type", 1);		//dilation
		s_smoothProgram->Enable();
		for (int l = 0; l < res_l; l++)
		{
			s_smoothProgram->SetFragmentParam1f("l", l+1);
			gridLayers[l]->DrawQuad();
		}
		renderedDensity->endRender();
		s_smoothProgram->Disable();
		s_RenderQuadVertexProgram->Disable();
	}

	//compute density gradient and store in density.xyz
	s_RenderQuadVertexProgram->Enable();
	renderedDensity->startRender();
	s_denGraProgram->SetFragmentTexParam("densityMap", renderedDensity->getTextureID());
	//s_denGraProgram->SetFragmentTexParam("densityMap", texlist[2]);
	s_denGraProgram->Enable();
	for(int l = 0; l < res_l; l++)
		gridLayers[l]->DrawQuad();
	renderedDensity->endRender();
	s_denGraProgram->Disable();
	s_RenderQuadVertexProgram->Disable();

	//compute 2nd derivative of density
	hessian_r1->startRender();
	hessian_r2->addtarget(1);
	hessian_r3->addtarget(2);
	GLenum buffers[] = { GL_COLOR_ATTACHMENT0_EXT, GL_COLOR_ATTACHMENT1_EXT, GL_COLOR_ATTACHMENT2_EXT};
	glDrawBuffers(3, buffers);
	s_RenderQuadVertexProgram->Enable();
	s_denHessianProgram->SetFragmentTexParam("gradientMap", renderedDensity->getTextureID());
	s_denHessianProgram->Enable();
	for(int l = 0; l < res_l; l++)
		gridLayers[l]->DrawQuad();
	s_denHessianProgram->Disable();
	s_RenderQuadVertexProgram->Disable();
	hessian_r1->endRender();
}

void FluidGridSimulator::createPuffBillBoard()
{
	billBoard = new RenderTexture(200, 200, true, GL_LINEAR, GL_TEXTURE_2D);

	billBoard->startRender();

	glViewport(0, 0, billBoard->getResX(), billBoard->getResY());
	//glViewport(0, 0, vpW, vpH);
	glClearColor(0, 0, 0, 0);
	glClearDepth(1.f);
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(45, (float)billBoard->getResX() / (float)billBoard->getResY(), 0.01f, 30);
	//gluPerspective(45, (float)vpW / (float)vpH, 0.01f, 30);

	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	gluLookAt(0, 0, 10.f, 0, 0, 0, 0, 1, 0);


	glEnable(GL_DEPTH_TEST);glDepthMask(GL_TRUE);
	s_billboardProgram->Enable();
	glutSolidSphere(3, 500, 500);
	s_billboardProgram->Disable();
	glDisable(GL_DEPTH_TEST);glDepthMask(GL_FALSE);

	billBoard->endRender();
}

// rendering function
void FluidGridSimulator::render(int w, int h)
{	
	IRenderTarget* textureToDisplay;
	vpW = w;
	vpH = h;

	static int stime = 0;
	//if (stime > 150)
	//	setPause(true); 
/*
	if the animation is not stopped, smoothing the density data
*/
	if (!getPause())
	{
		smoothingSmokeData();
		/////
		stime++;
	}
	//else printf("time: %d\n", stime);

	/*
	if the window is resized, dimensions of textures (intersection, footprint) should be changed to have the same size
*/
	glClearColor(0, 0, 0, 0);
	if ((intersection == NULL) || (intersection->getResX() != w) || (intersection->getResY() != h))
	{
		delete intersection;
		delete footprint;
		delete marktex;
		delete initpointtex;
		delete tmpadvecttex;
		delete ndotetex;
		delete edgetex;
		delete depthtex;
		delete pointCollection;
		delete advectCollection;
		delete advectCollection_2v;
		delete pointCollection_2v;
		delete tmpadvecttex_2v;
		glDeleteLists(dispListId, 1);

		if (intersection != NULL)
		{
			delete[] randIntensity;
		}
		randIntensity = new GLfloat[w*h];
		for (int i = 0; i < w; i++)
		{
			for (int j = 0; j < h; j++)
			{
				randIntensity[i*h + j] = (float)rand()/RAND_MAX;
			}
		}
		glDeleteTextures(1, &randTex);
		glGenTextures(1, &randTex);
		glBindTexture(GL_TEXTURE_RECTANGLE, randTex);
		glTexParameterf(GL_TEXTURE_RECTANGLE, GL_TEXTURE_WRAP_S, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_RECTANGLE, GL_TEXTURE_WRAP_T, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_RECTANGLE, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
		glTexParameterf(GL_TEXTURE_RECTANGLE, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
		glTexImage2D(GL_TEXTURE_RECTANGLE, 0, GL_RGBA, w,
			h, 0, GL_RED, GL_FLOAT, randIntensity);
		glBindTexture(GL_TEXTURE_RECTANGLE, 0);

		genSamplingList();

		//initialize a display list
		dispListId = glGenLists(1);
		glNewList(dispListId, GL_COMPILE);
		glBegin(GL_POINTS);
		for (int x = 0; x < w; x++)
			for (int y = 0; y < h; y++)
				glVertex3f((float)x + 0.f, (float)y + 0.f, randIntensity[x*h + y]);
		glEnd();
		glEndList();

		intersection = new RenderTexture(w, h, false, GL_NEAREST);
		footprint = new RenderTexture(w, h, false, GL_NEAREST);
		renderedCurvature = new RenderTexture(w, h, GL_NEAREST);
		marktex = new RenderTexture(w, h, true, GL_NEAREST);
		ndotetex = new RenderTexture(w, h, false, GL_NEAREST);
		edgetex = new IteratedRenderTexture(w, h, false, GL_NEAREST);
		depthtex = new RenderTexture(w, h, false, GL_NEAREST, GL_TEXTURE_2D);
		pointCollection = new RenderTexture(w, h, false, GL_NEAREST);
		advectCollection = new RenderTexture(w, h, false, GL_NEAREST);
		initpointtex = new RenderTexture(w, h, true, GL_NEAREST);
		tmpadvecttex = new RenderTexture(w, h, true, GL_NEAREST);
		advectCollection_2v = new RenderTexture(w, h, false, GL_NEAREST);
		pointCollection_2v = new RenderTexture(w, h, false, GL_NEAREST);
		tmpadvecttex_2v = new RenderTexture(w, h, true, GL_NEAREST);
	}

	switch(displaymode)
	{
	case 0:
		//render scene to img buff
		renderFlat3DSliceQuad(w, h, renderedDensity, is2Dsmoothed);
		renderFlat3D(w, h);
		//createPuffBillBoard();
		return;
	case 1:
		renderFlat3DSliceQuad(w, h, renderedDensity, is2Dsmoothed);
		renderRealSmoke(w, h);
		renderFlat3DSliceQuad(w, h, density, 0.f);					//not smoothed
		renderRealSmoke(w, h, w, 0);
		return;
	case 2:
		renderFlat3DSliceQuad(w, h, renderedDensity, is2Dsmoothed);
		ContourRendering();
		renderTextures(w, h, edgetex);
		return;
	case 3:
		renderTextures(w, h, ndotetex);
		return;
	case 4:
		//renderRealSmoke(w, h, density);
		return;
	case 5:
		textureToDisplay = density;
		break;
	case 6:
		textureToDisplay = vorticity;
		break;
	case 7:
		textureToDisplay = vorticityForce;
		break;
	case 8:
		textureToDisplay = temperature;
		break;
	case 9:
		renderGridPoint(w, h);
		return;
	}

	s_RenderQuadVertexProgram->Enable();

	renderTexturedQuad(textureToDisplay, 0.5, 0.5f);

	s_RenderQuadVertexProgram->Disable();
}

// calculate the size of the quads used during slicing 3D volume rendering
void calculateQuadSize(GLdouble ModelViewMatrix[16], float size[3], float pos[3], float* resultSize, float* resultCenter)
{

	float vert1[] = {pos[0]-size[0], pos[1]-size[1], pos[2]-size[2]};
	float vert2[] = {pos[0]-size[0], pos[1]-size[1], pos[2]+size[2]};
	float vert3[] = {pos[0]-size[0], pos[1]+size[1], pos[2]-size[2]};
	float vert4[] = {pos[0]-size[0], pos[1]+size[1], pos[2]+size[2]};
	float vert5[] = {pos[0]+size[0], pos[1]-size[1], pos[2]-size[2]};
	float vert6[] = {pos[0]+size[0], pos[1]-size[1], pos[2]+size[2]};
	float vert7[] = {pos[0]+size[0], pos[1]+size[1], pos[2]-size[2]};
	float vert8[] = {pos[0]+size[0], pos[1]+size[1], pos[2]+size[2]};

	float* vertices[] = {vert1, vert2, vert3, vert4, vert5, vert6, vert7, vert8};

	GLdouble rpos[4];
	float minX, maxX, minY, maxY, minZ, maxZ;
	GLdouble ProjMatrix[]={1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1};
	GLint viewport[] = {0,0,1,1};
	gluProject(vert1[0],
		vert1[1],
		vert1[2],
		ModelViewMatrix,
		ProjMatrix,
		viewport,
		&rpos[0],
		&rpos[1],
		&rpos[2]);


	minX = maxX = rpos[0]; minY = maxY = rpos[1]; minZ = maxZ = rpos[2];

	for(int i = 1; i < 8; i++)
	{
		gluProject(vertices[i][0],
			vertices[i][1],
			vertices[i][2],
			ModelViewMatrix,
			ProjMatrix,
			viewport,
			&rpos[0],
			&rpos[1],
			&rpos[2]);
		if(rpos[0] < minX) minX = rpos[0];
		if(rpos[0] > maxX) maxX = rpos[0];
		if(rpos[1] < minY) minY = rpos[1];
		if(rpos[1] > maxY) maxY = rpos[1];
		if(rpos[2] < minZ) minZ = rpos[2];
		if(rpos[2] > maxZ) maxZ = rpos[2];
	}
	maxX = maxX * 2.0 - 1.0;
	maxY = maxY * 2.0 - 1.0;
	maxZ = maxZ * 2.0 - 1.0;
	minX = minX * 2.0 - 1.0;
	minY = minY * 2.0 - 1.0;
	minZ = minZ * 2.0 - 1.0;

	resultSize[0] = (maxX - minX) / 2.0;
	resultSize[1] = (maxY - minY) / 2.0;
	resultSize[2] = (maxZ - minZ) / 2.0;
	resultCenter[0] = (maxX + minX) / 2.0;
	resultCenter[1] = (maxY + minY) / 2.0;
	resultCenter[2] = (maxZ + minZ) / 2.0;
}

void FluidGridSimulator::calculateMatrix(float* light, GLdouble* lmvp, GLdouble* emvp, GLdouble* hwmv, bool &f2b)
{
	//model view proj matrix of the eye
	glGetDoublev(GL_MODELVIEW_MATRIX, emvp);
	glMatrixMode(GL_PROJECTION); glPushMatrix();
	glMultMatrixd(emvp);
	glGetDoublev(GL_PROJECTION_MATRIX, emvp);
	glMatrixMode(GL_PROJECTION); glPopMatrix();

	//model view proj matrix of the light
	glViewport(0, 0, lw, lh);
	glMatrixMode(GL_MODELVIEW); glPushMatrix();
	glLoadIdentity();
	//gluPerspective(90, lw / lh, 0.01f, 30.f);
	glOrtho(-5, 5, -5, 5, 0.01f, 30);
	gluLookAt(light[0], light[1], light[2], 0, 0, 0, 0, 1, 0);
	glGetDoublev(GL_MODELVIEW_MATRIX, lmvp);
	glMatrixMode(GL_MODELVIEW); glPopMatrix();

	//model view matrix of the halfway coordinate (currently, it's same as mv matrix of the eye)
	Matrix4f emv;
	glGetFloatv(GL_MODELVIEW_MATRIX, emv.data);				//model view matrix of the eye
	Vector3f edir(-emv.data[2], -emv.data[6], -emv.data[10]); edir.normalize();
	Vector3f ldir(-light[0], -light[1], -light[2]); ldir.normalize();
	float EdotL = edir.dotProduct(ldir);
	Vector3f bisector;
	if (EdotL >= 0)
	{
		bisector = edir + ldir;
		f2b = true;
	}
	else
	{
		bisector = -edir + ldir;
		f2b = false;
	}
	glMatrixMode(GL_MODELVIEW); glPushMatrix();
	glLoadIdentity();
	gluLookAt(0, 0, 0, bisector.x, bisector.y, bisector.z, 0, 1, 0);
	glGetDoublev(GL_MODELVIEW_MATRIX, hwmv);
	glMatrixMode(GL_MODELVIEW); glPopMatrix();

}

void FluidGridSimulator::renderBoundingBox(float hx, float hy, float hz)
{
	//x = -hx;
	glBegin(GL_QUADS);
	glVertex3f(-hx, -hy, -hz); // The bottom left corner
	glVertex3f(-hx, -hy, hz); // The bottom right corner
	glVertex3f(-hx, hy, hz); // The top right corner
	glVertex3f(-hx, hy, -hz); // The top left corner
	glEnd();
	// y = -hy;
	glBegin(GL_QUADS);
	glVertex3f(-hx, -hy, -hz); // The bottom left corner
	glVertex3f(hx, -hy, -hz); // The bottom right corner
	glVertex3f(hx, -hy, hz); // The top right corner
	glVertex3f(-hx, -hy, hz); // The top left corner
	glEnd();
	// z = -hz;
	glBegin(GL_QUADS);
	glVertex3f(-hx, -hy, -hz); // The bottom left corner
	glVertex3f(-hx, hy, -hz); // The bottom right corner
	glVertex3f(hx, hy, -hz); // The top right corner
	glVertex3f(hx, -hy, -hz); // The top left corner
	glEnd();
	// x = hx;
	glBegin(GL_QUADS);
	glVertex3f(hx, hy, hz); // The bottom left corner
	glVertex3f(hx, -hy, hz); // The bottom right corner
	glVertex3f(hx, -hy, -hz); // The top right corner
	glVertex3f(hx, hy, -hz); // The top left corner
	glEnd();
	// y = hy;
	glBegin(GL_QUADS);
	glVertex3f(-hx, hy, -hz); // The bottom left corner
	glVertex3f(-hx, hy, hz); // The bottom right corner
	glVertex3f(hx, hy, hz); // The top right corner
	glVertex3f(hx, hy, -hz); // The top left corner
	glEnd();
	// z = hz;
	glBegin(GL_QUADS);
	glVertex3f(-hx, hy, hz); // The bottom left corner
	glVertex3f(-hx, -hy, hz); // The bottom right corner
	glVertex3f(hx, -hy, hz); // The top right corner
	glVertex3f(hx, hy, hz); // The top left corner
	glEnd();
}


inline void FullScreenQuad()
{
	glMatrixMode(GL_MODELVIEW);
	glPushMatrix();
	glLoadIdentity();
	glMatrixMode(GL_PROJECTION);
	glPushMatrix();
	glLoadIdentity();
	glBegin(GL_QUADS);
	glVertex3f(-1,-1,0);  glVertex3f(1,-1,0);  glVertex3f(1,1,0);  glVertex3f(-1,1,0);
	glEnd();
	glPopMatrix();
	glMatrixMode(GL_MODELVIEW);
	glPopMatrix();
}

// render density in 3D using slicing quad
void FluidGridSimulator::renderFlat3DSliceQuad(int w, int h, IRenderTarget* volumedata, int isSmoothed)
{
	///////////
	// set GPU program attributes, compute slicing quad position and size, render quads
	/////////

	GLdouble lightModelViewProj[16];
	GLdouble eyeModelViewProj[16];
	GLdouble halfwayModelView[16];
	//float light[] = {0.5, 0.3, -0.5f};//{-0.5, -0.3, 0.5f};//{0.3f, 0.5f, 2.f};
	bool f2b;
	calculateMatrix(lightPos, lightModelViewProj, eyeModelViewProj, halfwayModelView, f2b);

	float size[3];
	float center[3];
	float m_size[]={sizeX / 2, sizeZ / 2, sizeY / 2};
	float m_pos[]={posX,posZ,posY};
	calculateQuadSize(halfwayModelView, m_size, m_pos, size, center);

	//int numSlices = renderSliceCount;
	//if(renderSliceCount == 0)
	//	numSlices =  res_m * 1.5;
	int numSlices = res_m * 5;

	float z = size[2];
	float step = size[2] * 2 / (numSlices - 1);

	float miplevel = 1.f;
	int nw = w / miplevel;
	int nh = h / miplevel;

	if ((imgbuff == NULL) || (imgbuff->getResX() != nw) || (imgbuff->getResY() != nh))
	{
		delete imgbuff;
		imgbuff = new IteratedRenderTexture(nw, nh, false, GL_LINEAR, GL_TEXTURE_2D);
	}

	//initialize light buffer to light color
	lightbuffer->startRender();
	glClearColor(1, 1, 1, 0);
	glClear(GL_COLOR_BUFFER_BIT);
	lightbuffer->endRender();

	imgbuff->startRender();
	glViewport(0, 0, nw, nh);
	glClearColor(0, 0, 0, 0);
	glClear(GL_COLOR_BUFFER_BIT);
	imgbuff->endRender();

	//set model view matrix to halfway coordinate
	glMatrixMode(GL_MODELVIEW); glPushMatrix();
	glLoadIdentity();
	glMultMatrixd(halfwayModelView);

	glEnable(GL_BLEND);
	glBlendEquation(GL_FUNC_ADD);
	GLenum sfactor = GL_ONE_MINUS_DST_ALPHA;
	GLenum dfactor = GL_ONE;
	if (!f2b)
	{
		sfactor = GL_ONE;
		dfactor = GL_ONE_MINUS_SRC_ALPHA;
	}
	for(int i = 0; i < numSlices; i++)
	{
		//pass 1
		imgbuff->startRender();
		glViewport(0, 0, nw, nh);
		glBlendFunc(sfactor, dfactor);
		s_VolumeRenderSliceFlat3DProgram->SetVertexParam3f("gridscale", 1.0 / m_size[0], 1.0 / m_size[1], 1.0 / m_size[2]);
		s_VolumeRenderSliceFlat3DProgram->SetVertexParam3f("gridcenter", m_pos[0], m_pos[1], m_pos[2]);
		s_VolumeRenderSliceFlat3DProgram->SetVertexMatrixParam("ModelViewProjCoord", eyeModelViewProj);
		s_VolumeRenderSliceFlat3DProgram->SetFragmentParam3f("resolution", res_m, res_n, res_l);
		s_VolumeRenderSliceFlat3DProgram->SetFragmentParam2f("lvpsize", lw, lh);
		s_VolumeRenderSliceFlat3DProgram->SetFragmentTexParam("density", volumedata->getTextureID());
		s_VolumeRenderSliceFlat3DProgram->SetFragmentTexParam("offsetTexture", layerOffsetTexture);
		s_VolumeRenderSliceFlat3DProgram->SetFragmentMatrixParam("lightModelViewProj", lightModelViewProj);
		s_VolumeRenderSliceFlat3DProgram->SetFragmentTexParam("lightBuffer", lightbuffer->getTextureID());
		//s_VolumeRenderSliceFlat3DProgram->SetFragmentTexParam("positionMap", footprint->getTextureID());
		s_VolumeRenderSliceFlat3DProgram->Enable();
		glBegin(GL_QUADS);
		glTexCoord3f(-size[0],-size[1],z);
		glVertex3fv(center);
		glTexCoord3f(size[0],-size[1],z);
		glVertex3fv(center);
		glTexCoord3f(size[0],size[1],z);
		glVertex3fv(center);
		glTexCoord3f(-size[0],size[1],z);
		glVertex3fv(center);
		glEnd();
		s_VolumeRenderSliceFlat3DProgram->Disable();	
		imgbuff->endRender();

		//pass 2
		lightbuffer->startRender();
		glBlendFunc(GL_ONE, GL_ONE_MINUS_SRC_ALPHA);
		s_lightProgram->SetVertexParam3f("gridscale", 1.0 / m_size[0], 1.0 / m_size[1], 1.0 / m_size[2]);
		s_lightProgram->SetVertexParam3f("gridcenter", m_pos[0], m_pos[1], m_pos[2]);
		s_lightProgram->SetVertexMatrixParam("ModelViewProjCoord", lightModelViewProj);
		s_lightProgram->SetFragmentParam3f("resolution", res_m, res_n, res_l);
		s_lightProgram->SetFragmentTexParam("density", volumedata->getTextureID());
		s_lightProgram->SetFragmentTexParam("offsetTexture", layerOffsetTexture);
		s_lightProgram->Enable();
		glBegin(GL_QUADS);
		glTexCoord3f(-size[0],-size[1],z);
		glVertex3fv(center);
		glTexCoord3f(size[0],-size[1],z);
		glVertex3fv(center);
		glTexCoord3f(size[0],size[1],z);
		glVertex3fv(center);
		glTexCoord3f(-size[0],size[1],z);
		glVertex3fv(center);
		glEnd();
		s_lightProgram->Disable();
		lightbuffer->endRender();

		z -= step;
	}
	glMatrixMode(GL_MODELVIEW); glPopMatrix();
	glDisable(GL_BLEND);

	//detect smoke surface
	SurfaceDetect();

	if (isSmoothed)
	{
		//smooth imgbuff
		imgbuff->startRender();
		glViewport(0, 0, w, h);
		s_RenderQuadVertexProgram->Enable();
		s_smooth2DProgram->SetFragmentTexParam("image", imgbuff->getTextureID());
		s_smooth2DProgram->SetFragmentTexParam("positionMap", footprint->getTextureID());
		s_smooth2DProgram->SetFragmentParam2f("scrSize", vpW, vpH);
		s_smooth2DProgram->Enable();
		glBegin(GL_QUADS);
		glTexCoord2f(0,0); glVertex3f(-1,-1,0);
		glTexCoord2f(1,0); glVertex3f(1,-1,0);
		glTexCoord2f(1,1); glVertex3f(1,1,0);
		glTexCoord2f(0,1); glVertex3f(-1,1,0);		
		glEnd();
		s_smooth2DProgram->Disable();
		s_RenderQuadVertexProgram->Disable();
		imgbuff->endRender();
	}

	//render background+smoke to imgbuff
	imgbuff->startRender();
	glViewport(0, 0, w, h);
	s_RenderQuadVertexProgram->Enable();
	s_realRenderProgram->SetFragmentTexParam("colorbuff", imgbuff->getTextureID());
	s_realRenderProgram->SetFragmentParam4f("bkgcolor", .35f,.35f,.35f,0.f);//0.5f, 0.5f, 0.5f, 0);
	s_realRenderProgram->SetFragmentParam1f("intensitymode", intensitymode);
	s_realRenderProgram->SetFragmentParam1f("a", gausssig);
	s_realRenderProgram->Enable();
	s_realRenderProgram->EnableFragmentTexParam("colorbuff");
	glBegin(GL_QUADS);
	glTexCoord2f(0,0); glVertex3f(-1,-1,0);
	glTexCoord2f(1,0); glVertex3f(1,-1,0);
	glTexCoord2f(1,1); glVertex3f(1,1,0);
	glTexCoord2f(0,1); glVertex3f(-1,1,0);		
	glEnd();
	s_realRenderProgram->DisableFragmentTexParam("colorbuff");
	s_realRenderProgram->Disable();
	s_RenderQuadVertexProgram->Disable();
	imgbuff->endRender();
}

void FluidGridSimulator::renderRealSmoke(int w, int h, int w0, int h0)
{
	//render to screen with img buff as a texture
	s_RenderQuadVertexProgram->Enable();
	s_realRenderProgram->SetFragmentTexParam("colorbuff", imgbuff->getTextureID());
	s_realRenderProgram->SetFragmentParam4f("bkgcolor", 0, 0, 0.f, 0);
	s_realRenderProgram->SetFragmentParam1f("intensitymode", 0.f);
	s_realRenderProgram->SetFragmentTexParam("ndoteMap", ndotetex->getTextureID());
	s_realRenderProgram->SetFragmentTexParam("toneMap", halftoneBand);
	s_realRenderProgram->Enable();
	s_realRenderProgram->EnableFragmentTexParam("colorbuff");
	glViewport(w0, h0, w, h);
	glBegin(GL_QUADS);
	glTexCoord2f(0,0); glVertex3f(-1,-1,0);
	glTexCoord2f(1,0); glVertex3f(1,-1,0);
	glTexCoord2f(1,1); glVertex3f(1,1,0);
	glTexCoord2f(0,1); glVertex3f(-1,1,0);		
	glEnd();
	s_realRenderProgram->DisableFragmentTexParam("colorbuff");
	s_realRenderProgram->Disable();
	s_RenderQuadVertexProgram->Disable();

	//render bounding cube
	glEnable(GL_BLEND);
	glBlendFunc(GL_ONE, GL_ZERO);
	if(renderBoundingCube)
	{	
		glColor4f(0,0,0,0);
		glPushMatrix();
		glTranslatef(posX, posZ, posY);
		glScalef(sizeX, sizeZ, sizeY);
		glutWireCube(1);
		glPopMatrix();
	}
	glDisable(GL_BLEND);
}

void FluidGridSimulator::renderTextures(int w, int h, IRenderTarget* tex)
{
	//render to screen with img buff as a texture
	s_RenderQuadVertexProgram->Enable();
	s_renderTexProgram->SetFragmentTexParam("colorbuff", tex->getTextureID());
	s_renderTexProgram->SetFragmentParam4f("bkgcolor", 0, 0, 1, 0);
	s_renderTexProgram->Enable();
	s_renderTexProgram->EnableFragmentTexParam("colorbuff");
	glViewport(0, 0, w, h);
	glBegin(GL_QUADS);
	glVertex3f(-1,-1,0);
	glVertex3f(1,-1,0);
	glVertex3f(1,1,0);
	glVertex3f(-1,1,0);		
	glEnd();
	s_renderTexProgram->DisableFragmentTexParam("colorbuff");
	s_renderTexProgram->Disable();
	s_RenderQuadVertexProgram->Disable();

	//render bounding cube
	glEnable(GL_BLEND);
	glBlendFunc(GL_ONE, GL_ZERO);
	if(renderBoundingCube)
	{	
		glColor4f(0,0,0,0);
		glPushMatrix();
		glTranslatef(posX, posZ, posY);
		glScalef(sizeX, sizeZ, sizeY);
		glutWireCube(1);
		glPopMatrix();
	}
	glDisable(GL_BLEND);
}

void FluidGridSimulator::SurfaceDetect()
{
		intersection->startRender();
			//render the back face
			s_IntersecbackProgram->Enable();
				glClear(GL_COLOR_BUFFER_BIT);
				glEnable(GL_CULL_FACE);
				glCullFace(GL_FRONT);
				glPushMatrix();
				glTranslatef(posX, posZ, posY);
				renderBoundingBox(sizeX / 2.f, sizeZ / 2.f, sizeY / 2.f);
				glPopMatrix();
				glDisable(GL_CULL_FACE);
			s_IntersecbackProgram->Disable();

			//render the front face, applying blending
			s_IntersecfrontProgram->Enable();
				glEnable(GL_BLEND);
				glBlendFunc(GL_ONE, GL_ONE);
				glBlendEquation(GL_FUNC_SUBTRACT);
				glEnable(GL_CULL_FACE);
				glCullFace(GL_BACK);
				glPushMatrix();
				glTranslatef(posX, posZ, posY);
				renderBoundingBox(sizeX / 2.f, sizeZ / 2.f, sizeY / 2.f);
				glPopMatrix();		
				glDisable(GL_CULL_FACE);
				glDisable(GL_BLEND);
			s_IntersecfrontProgram->Disable();
		intersection->endRender();

	/*
		locate the ray-smoke intersection coordinates (in object space) and store it to footprint texture
	*/
		ndotetex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f);
		glClear(GL_COLOR_BUFFER_BIT);
		ndotetex->endRender();

		footprint->startRender();
		ndotetex->addtarget(1);
		depthtex->addtarget(2);
		GLenum buffers[] = { GL_COLOR_ATTACHMENT0_EXT, GL_COLOR_ATTACHMENT1_EXT, GL_COLOR_ATTACHMENT2_EXT};
		glDrawBuffers(3, buffers);
		glClearColor(0.f, 0.f, 0.f, 0.f);
		glClear(GL_COLOR_BUFFER_BIT);
		s_RenderQuadVertexProgram->Enable();
		s_LocateIntersectProgram->SetFragmentTexParam("intersecMap", intersection->getTextureID());
		s_LocateIntersectProgram->SetFragmentTexParam("densityMap", renderedDensity->getTextureID());
		s_LocateIntersectProgram->SetFragmentTexParam("offsetTexture", layerOffsetTexture);
		s_LocateIntersectProgram->SetFragmentParam1f("cellsize", sizeX / res_m);
		s_LocateIntersectProgram->SetFragmentParam3f("gridscale", 1.0 / m_size[0], 1.0 / m_size[1], 1.0 / m_size[2]);
		s_LocateIntersectProgram->SetFragmentParam3f("gridcenter", m_pos[0], m_pos[1], m_pos[2]);
		s_LocateIntersectProgram->SetFragmentParam3f("resolution", res_m, res_n, res_l);
		s_LocateIntersectProgram->SetFragmentTexParam("imgbuff", imgbuff->getTextureID());
		s_LocateIntersectProgram->SetFragmentParam2f("scrSize", vpW, vpH);
		s_LocateIntersectProgram->SetFragmentParam1f("detectmode", detectMode);
		s_LocateIntersectProgram->SetFragmentParam1f("DENSITY_THRES", densityThres);
		s_LocateIntersectProgram->SetFragmentParam1f("INTENSITY_THRES", intensityThres);

		//s_LocateIntersectProgram->SetFragmentTexParam("hesrow1", hessian_r1->getTextureID());
		//s_LocateIntersectProgram->SetFragmentTexParam("hesrow2", hessian_r2->getTextureID());
		//s_LocateIntersectProgram->SetFragmentTexParam("hesrow3", hessian_r3->getTextureID());

		s_LocateIntersectProgram->Enable();

			glEnable(GL_BLEND);
			glBlendFunc(GL_ONE, GL_ONE);
			glBlendEquation(GL_FUNC_ADD);
			glViewport(0, 0, vpW, vpH);
			glBegin(GL_QUADS);
			glTexCoord2f(0,0); glVertex3f(-1,-1,0);
			glTexCoord2f(vpW,0); glVertex3f(1,-1,0);
			glTexCoord2f(vpW,vpH); glVertex3f(1,1,0);
			glTexCoord2f(0,vpH); glVertex3f(-1,1,0);		
			glEnd();
			glDisable(GL_BLEND);

		s_LocateIntersectProgram->Disable();
		s_RenderQuadVertexProgram->Disable();
		footprint->endRender();
}

/* 
At each time step, after compute simulation data, this function is called for rendering NPR smoke
*/
void FluidGridSimulator::renderFlat3D(int w, int h) 
{
/*
	if the camera was not changed and animation is paused, just re-render previous frame
*/
	if (getPause() || isCamChanged)
	{
		//RorateSmoke();
		//return;
	}

/*
	render back face and front face of the bounding box to intersection texture
	- By rendering front face, we can compute the ray-box intersection and eye-front face distance
	- By rendering back face while applying subtract blending, we can compute how far the ray travel 
	through the bounding box
*/
	switch (nprStyle)
	{
	case 0:
		NPRcoloring();
		break;
	case 1:
		Stippling();
		//SolidContour();
		return;
	case 2:
		Hatching();
		break;
	}

/*
	render bounding box
*/
	glEnable(GL_BLEND);
	glBlendFunc(GL_ONE, GL_ZERO);
	if(renderBoundingCube)
	{	
		glColor4f(0.8,0.8,0.8,0);
		glPushMatrix();
		glTranslatef(posX, posZ, posY);
		glScalef(sizeX, sizeZ, sizeY);
		glutWireCube(1);
		glPopMatrix();
		/////
		glColor4f(1.f, 0.f, 0.f, 0.f);
		glPushMatrix();
		glTranslatef(0.f, (m_windsource[2] - res_l*0.5f ) / res_m, (m_windsource[1] - res_n*0.5f) / res_m);
		glutWireCube(0.1f);
		glPopMatrix();	
	}
	glDisable(GL_BLEND);
}

void FluidGridSimulator::RorateSmoke()
{
	{
		isCamChanged = false;

		intersection->startRender();
			//render the back face
			s_IntersecbackProgram->Enable();
				glClear(GL_COLOR_BUFFER_BIT);
				glEnable(GL_CULL_FACE);
				glCullFace(GL_FRONT);
				glPushMatrix();
				glTranslatef(posX, posZ, posY);
				renderBoundingBox(sizeX / 2.f, sizeZ / 2.f, sizeY / 2.f);
				glPopMatrix();
				glDisable(GL_CULL_FACE);
			s_IntersecbackProgram->Disable();

			//render the front face, applying blending
			s_IntersecfrontProgram->Enable();
				glEnable(GL_BLEND);
				glBlendFunc(GL_ONE, GL_ONE);
				glBlendEquation(GL_FUNC_SUBTRACT);
				glEnable(GL_CULL_FACE);
				glCullFace(GL_BACK);
				glPushMatrix();
				glTranslatef(posX, posZ, posY);
				renderBoundingBox(sizeX / 2.f, sizeZ / 2.f, sizeY / 2.f);
				glPopMatrix();		
				glDisable(GL_CULL_FACE);
				glDisable(GL_BLEND);
			s_IntersecfrontProgram->Disable();
		intersection->endRender();

	/*
		locate the ray-smoke intersection coordinates (in object space) and store it to footprint texture
	*/
		ndotetex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f);
		glClear(GL_COLOR_BUFFER_BIT);
		ndotetex->endRender();

		footprint->startRender();
		ndotetex->addtarget(1);
		depthtex->addtarget(2);
		GLenum buffers[] = { GL_COLOR_ATTACHMENT0_EXT, GL_COLOR_ATTACHMENT1_EXT, GL_COLOR_ATTACHMENT2_EXT};
		glDrawBuffers(3, buffers);
		glClearColor(0.f, 0.f, 0.f, 0.f);
		glClear(GL_COLOR_BUFFER_BIT);
		s_RenderQuadVertexProgram->Enable();
		s_LocateIntersectProgram->SetFragmentTexParam("intersecMap", intersection->getTextureID());
		s_LocateIntersectProgram->SetFragmentTexParam("densityMap", renderedDensity->getTextureID());
		s_LocateIntersectProgram->SetFragmentTexParam("offsetTexture", layerOffsetTexture);
		s_LocateIntersectProgram->SetFragmentParam1f("cellsize", sizeX / res_m);
		s_LocateIntersectProgram->SetFragmentParam3f("gridscale", 1.0 / m_size[0], 1.0 / m_size[1], 1.0 / m_size[2]);
		s_LocateIntersectProgram->SetFragmentParam3f("gridcenter", m_pos[0], m_pos[1], m_pos[2]);
		s_LocateIntersectProgram->SetFragmentParam3f("resolution", res_m, res_n, res_l);
		s_LocateIntersectProgram->Enable();

			glEnable(GL_BLEND);
			glBlendFunc(GL_ONE, GL_ONE);
			glBlendEquation(GL_FUNC_ADD);
			glViewport(0, 0, vpW, vpH);
			glBegin(GL_QUADS);
			glTexCoord2f(0,0); glVertex3f(-1,-1,0);
			glTexCoord2f(vpW,0); glVertex3f(1,-1,0);
			glTexCoord2f(vpW,vpH); glVertex3f(1,1,0);
			glTexCoord2f(0,vpH); glVertex3f(-1,1,0);		
			glEnd();
			glDisable(GL_BLEND);

		s_LocateIntersectProgram->Disable();
		s_RenderQuadVertexProgram->Disable();
		footprint->endRender();
	}

	//render scene to img buff
	renderFlat3DSliceQuad(vpW, vpH, renderedDensity);

/*
	//stippling
	//make an initial point distribution and mark
	{
		initpointtex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_intertexturesProgram->SetVertexTexParam("positionMap", footprint->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("imgbuff", imgbuff->getTextureID());
		s_intertexturesProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_intertexturesProgram->Enable();

		glCallList(samplingListId);

		s_intertexturesProgram->Disable();
		initpointtex->endRender();
	}

	//mark covering region
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		marktex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glClearDepth(1.f);glClear(GL_DEPTH_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_markRegionProgram->SetVertexTexParam("initpointMap", initpointtex->getTextureID());
		s_markRegionProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_markRegionProgram->Enable();

		glCallList(samplingListId);					

		s_markRegionProgram->Disable();
		marktex->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}

	//advect point collection
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		initpointtex->startRender();
		glViewport(0, 0, vpW, vpH);
		glClearDepth(1.f);glClear(GL_DEPTH_BUFFER_BIT);
		s_rotateSampPoint->SetVertexTexParam("pointMap", pointCollection->getTextureID());
		s_rotateSampPoint->SetVertexTexParam("markMap", marktex->getTextureID());
		s_rotateSampPoint->SetVertexParam2f("scrSize", vpW, vpH);
		s_rotateSampPoint->Enable();
		glCallList(dispListId);
		s_rotateSampPoint->Disable();
		initpointtex->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}

	//distribute points
	{
		pointCollection->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_distPointProgram->SetVertexTexParam("initpointMap", initpointtex->getTextureID());
		s_distPointProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_distPointProgram->SetFragmentTexParam("positionMap", footprint->getTextureID());
		s_distPointProgram->Enable();

		glCallList(samplingListId);

		s_distPointProgram->Disable();
		pointCollection->endRender();
	}

	//hatching with by accessing the texture stored curvature info
	{
		glClearColor(1.f, 1.f, 1.f, 1.f); glClear(GL_COLOR_BUFFER_BIT);
		glEnable(GL_POINT_SMOOTH);
		glPointSize(1.f);
		glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
		glEnable(GL_BLEND);
		glBlendFunc(GL_DST_COLOR, GL_ZERO);
		glViewport(0, 0, vpW, vpH);
		s_stippleProgram->SetVertexTexParam("pointMap", pointCollection->getTextureID());
		s_stippleProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_stippleProgram->Enable();

		glCallList(dispListId);

		s_stippleProgram->Disable();
		glDisable(GL_BLEND);
		glDisable(GL_POINT_SMOOTH);
		glPointSize(1.f);
	}
*/

	//edge detection filtering
	{
		edgetex->startRender();
		s_RenderQuadVertexProgram->Enable();
		s_edgedetectProgram->SetFragmentTexParam("ndotetext", ndotetex->getTextureID());
		s_edgedetectProgram->SetFragmentTexParam("imgbuff", imgbuff->getTextureID());
		s_edgedetectProgram->SetFragmentParam2f("scrSize", vpW, vpH);
		s_edgedetectProgram->SetFragmentParam1f("edgemode", 0);
		s_edgedetectProgram->Enable();
		glViewport(0, 0, vpW, vpH);
		glBegin(GL_QUADS);
		glTexCoord2f(0,0); glVertex3f(-1,-1,0);
		glTexCoord2f(1,0); glVertex3f(1,-1,0);
		glTexCoord2f(1,1); glVertex3f(1,1,0);
		glTexCoord2f(0,1); glVertex3f(-1,1,0);		
		glEnd();
		s_edgedetectProgram->Disable();
		s_RenderQuadVertexProgram->Disable();
		edgetex->endRender();
	}

	if (edgemode == 1)
		//combine edges
	{
		edgetex->startRender();
		s_RenderQuadVertexProgram->Enable();
		s_edgedetectProgram->SetFragmentTexParam("ndotetext", ndotetex->getTextureID());
		s_edgedetectProgram->SetFragmentTexParam("imgbuff", imgbuff->getTextureID());
		s_edgedetectProgram->SetFragmentParam2f("scrSize", vpW, vpH);
		s_edgedetectProgram->SetFragmentParam1f("edgemode", 1);
		s_edgedetectProgram->SetFragmentParam1f("joinmode", joinmode);
		s_edgedetectProgram->SetFragmentTexParam("edgeMap", edgetex->getTextureID());
		s_edgedetectProgram->SetFragmentParam1f("a", gausssig);
		s_edgedetectProgram->Enable();
		glViewport(0, 0, vpW, vpH);
		glBegin(GL_QUADS);
		glTexCoord2f(0,0); glVertex3f(-1,-1,0);
		glTexCoord2f(1,0); glVertex3f(1,-1,0);
		glTexCoord2f(1,1); glVertex3f(1,1,0);
		glTexCoord2f(0,1); glVertex3f(-1,1,0);		
		glEnd();
		s_edgedetectProgram->Disable();
		s_RenderQuadVertexProgram->Disable();
		edgetex->endRender();
	}

	////////////////////hatching
	//make an initial point distribution and mark
	{
		initpointtex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_intertexturesProgram->SetVertexTexParam("positionMap", footprint->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("imgbuff", imgbuff->getTextureID());
		s_intertexturesProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_intertexturesProgram->Enable();

		glCallList(samplingListId);					//need to be improved (i+=4; j+=4)

		s_intertexturesProgram->Disable();
		initpointtex->endRender();
	}

	//mark covering region
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		marktex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glClearDepth(1.f);glClear(GL_DEPTH_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_markRegionProgram->SetVertexTexParam("initpointMap", initpointtex->getTextureID());
		s_markRegionProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_markRegionProgram->Enable();

		glCallList(samplingListId);					

		s_markRegionProgram->Disable();
		marktex->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}

	//advect point collection
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		initpointtex->startRender();
		glViewport(0, 0, vpW, vpH);
		glClearDepth(1.f);glClear(GL_DEPTH_BUFFER_BIT);
		s_rotateSampPoint->SetVertexTexParam("pointMap", pointCollection->getTextureID());
		s_rotateSampPoint->SetVertexTexParam("markMap", marktex->getTextureID());
		s_rotateSampPoint->SetVertexParam2f("scrSize", vpW, vpH);
		s_rotateSampPoint->Enable();
		glCallList(dispListId);
		s_rotateSampPoint->Disable();
		initpointtex->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}

	//distribute points
	{
		pointCollection->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_distPointProgram->SetVertexTexParam("initpointMap", initpointtex->getTextureID());
		s_distPointProgram->Enable();

		glCallList(samplingListId);

		s_distPointProgram->Disable();
		pointCollection->endRender();
	}

	//hatching with by accessing the texture stored curvature info
	{
		glClearColor(1.f, 1.f, 1.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
		glEnable(GL_BLEND);
		glBlendFunc(GL_DST_COLOR, GL_ZERO);
		glEnable(GL_ALPHA_TEST);
		glAlphaFunc(GL_LESS, 0.8f);
		glViewport(0, 0, vpW, vpH);
		s_hatchingProgram->SetVertexTexParam("positionMap", footprint->getTextureID());
		s_hatchingProgram->SetVertexTexParam("pointMap", pointCollection->getTextureID());
		s_hatchingProgram->SetVertexTexParam("densityMap", renderedDensity->getTextureID());
		s_hatchingProgram->SetVertexTexParam("offsetTexture", layerOffsetTexture);
		s_hatchingProgram->SetVertexTexParam("hesrow1", hessian_r1->getTextureID());
		s_hatchingProgram->SetVertexTexParam("hesrow2", hessian_r2->getTextureID());
		s_hatchingProgram->SetVertexTexParam("hesrow3", hessian_r3->getTextureID());
		s_hatchingProgram->SetVertexParam3f("gridscale", 1.0 / m_size[0], 1.0 / m_size[1], 1.0 / m_size[2]);
		s_hatchingProgram->SetVertexParam3f("gridcenter", m_pos[0], m_pos[1], m_pos[2]);
		s_hatchingProgram->SetVertexParam3f("resolution", res_m, res_n, res_l);
		s_hatchingProgram->SetVertexParam1f("recsize", recsize);
		s_hatchingProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_hatchingProgram->SetVertexParam1f("hatchmode", hatchmode);
		s_hatchingProgram->SetVertexTexParam("edgeMap", edgetex->getTextureID());
		s_hatchingProgram->SetFragmentTexParam("strokeMap", strokeTexture);
		s_hatchingProgram->SetGeometryTexParam("randMap", randTex);
		s_hatchingProgram->Enable();

		glCallList(dispListId);

		s_hatchingProgram->Disable();
		glDisable(GL_BLEND);
		glDisable(GL_ALPHA_TEST);
	}

	//render bounding box
	glEnable(GL_BLEND);
	glBlendFunc(GL_ONE, GL_ZERO);
	if(renderBoundingCube)
	{	
		glColor4f(0.8,0.8,0.8,0);
		glPushMatrix();
		glTranslatef(posX, posZ, posY);
		glScalef(sizeX, sizeZ, sizeY);
		glutWireCube(1);
		glPopMatrix();
	}
	glDisable(GL_BLEND);

}

void FluidGridSimulator::ContourRendering()
{
	//edge detection filtering
	{
		edgetex->startRender();
		s_RenderQuadVertexProgram->Enable();
		s_edgedetectProgram->SetFragmentTexParam("ndotetext", ndotetex->getTextureID());
		s_edgedetectProgram->SetFragmentTexParam("imgbuff", imgbuff->getTextureID());
		s_edgedetectProgram->SetFragmentTexParam("randMap", randTex);
		s_edgedetectProgram->SetFragmentTexParam("edgeMap", edgetex->getTextureID());
		s_edgedetectProgram->SetFragmentParam2f("scrSize", vpW, vpH);
		s_edgedetectProgram->SetFragmentParam1f("edgemode", 0);
		s_edgedetectProgram->SetFragmentParam1f("a", gausssig);
		s_edgedetectProgram->Enable();
		glViewport(0, 0, vpW, vpH);
		glBegin(GL_QUADS);
		glTexCoord2f(0,0); glVertex3f(-1,-1,0);
		glTexCoord2f(1,0); glVertex3f(1,-1,0);
		glTexCoord2f(1,1); glVertex3f(1,1,0);
		glTexCoord2f(0,1); glVertex3f(-1,1,0);		
		glEnd();
		s_edgedetectProgram->Disable();
		s_RenderQuadVertexProgram->Disable();
		edgetex->endRender();
	}

	if (edgemode == 1)
	//combine edges
	{
		edgetex->startRender();
		s_RenderQuadVertexProgram->Enable();
		s_edgedetectProgram->SetFragmentTexParam("ndotetext", ndotetex->getTextureID());
		s_edgedetectProgram->SetFragmentTexParam("imgbuff", imgbuff->getTextureID());
		s_edgedetectProgram->SetFragmentParam2f("scrSize", vpW, vpH);
		s_edgedetectProgram->SetFragmentParam1f("edgemode", 1);
		s_edgedetectProgram->SetFragmentParam1f("joinmode", joinmode);
		s_edgedetectProgram->SetFragmentTexParam("edgeMap", edgetex->getTextureID());
		s_edgedetectProgram->SetFragmentParam1f("a", gausssig);
		s_edgedetectProgram->Enable();
		glViewport(0, 0, vpW, vpH);
		glBegin(GL_QUADS);
		glTexCoord2f(0,0); glVertex3f(-1,-1,0);
		glTexCoord2f(1,0); glVertex3f(1,-1,0);
		glTexCoord2f(1,1); glVertex3f(1,1,0);
		glTexCoord2f(0,1); glVertex3f(-1,1,0);		
		glEnd();
		s_edgedetectProgram->Disable();
		s_RenderQuadVertexProgram->Disable();
		edgetex->endRender();
	}

/*	///////////////////////////////////////////test
	//make an initial point distribution and mark
	{
		initpointtex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_intertexturesProgram->SetVertexTexParam("positionMap", footprint->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("imgbuff", imgbuff->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("ndoteMap", ndotetex->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("edgeMap", edgetex->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("toneMap", halftoneBand);
		s_intertexturesProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_intertexturesProgram->SetVertexParam1f("intensitymode", intensitymode);
		s_intertexturesProgram->SetVertexParam1f("a", gausssig);
		s_intertexturesProgram->Enable();

		//glCallList(dispListId);					//need to be improved (i+=4; j+=4)
		glCallList(samplingListId);

		s_intertexturesProgram->Disable();
		initpointtex->endRender();
	}

	//mark covering region
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		marktex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glClearDepth(1.f);glClear(GL_DEPTH_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_markRegionProgram->SetVertexTexParam("initpointMap", initpointtex->getTextureID());
		s_markRegionProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_markRegionProgram->Enable();

		glCallList(samplingListId);					

		s_markRegionProgram->Disable();
		marktex->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}

	//advect point collection
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		initpointtex->startRender();
		glViewport(0, 0, vpW, vpH);
		glClearDepth(1.f);glClear(GL_DEPTH_BUFFER_BIT);
		s_advectSampPoint->SetVertexTexParam("pointMap", pointCollection->getTextureID());
		s_advectSampPoint->SetVertexTexParam("positionMap", footprint->getTextureID());
		s_advectSampPoint->SetVertexTexParam("velocityMap", velocity->getTextureID());
		s_advectSampPoint->SetVertexTexParam("offsetTexture", layerOffsetTexture);
		s_advectSampPoint->SetVertexTexParam("markMap", marktex->getTextureID());
		s_advectSampPoint->SetVertexParam3f("gridscale", 1.0 / m_size[0], 1.0 / m_size[1], 1.0 / m_size[2]);
		s_advectSampPoint->SetVertexParam3f("gridcenter", m_pos[0], m_pos[1], m_pos[2]);
		s_advectSampPoint->SetVertexParam3f("resolution", res_m, res_n, res_l);
		s_advectSampPoint->SetVertexParam2f("scrSize", vpW, vpH);
		s_advectSampPoint->Enable();
		glCallList(dispListId);
		s_advectSampPoint->Disable();
		initpointtex->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}

	//distribute points
	{
		pointCollection->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_distPointProgram->SetVertexTexParam("initpointMap", initpointtex->getTextureID());
		s_distPointProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_distPointProgram->SetFragmentTexParam("positionMap", footprint->getTextureID());
		s_distPointProgram->Enable();

		glCallList(samplingListId);

		s_distPointProgram->Disable();
		pointCollection->endRender();
	}
	///////////////////////////////////////////////////////////////*/

	//texturing contours horizontally
	{
		glEnable(GL_ALPHA_TEST);
		glAlphaFunc(GL_LESS, 0.5f);
		glClearColor(1, 1, 1, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_texcontourProgram->SetGeometryTexParam("positionMap", footprint->getTextureID());
		s_texcontourProgram->SetGeometryTexParam("edgeMap", edgetex->getTextureID());
		s_texcontourProgram->SetFragmentTexParam("strokeMap", contourstrokeTexture);
		s_texcontourProgram->SetFragmentTexParam("depthMap", depthtex->getTextureID());
		s_texcontourProgram->SetGeometryParam2f("scrSize", vpW, vpH);
		s_texcontourProgram->SetGeometryParam1f("rsize", m_contourtexsize);
		s_texcontourProgram->SetGeometryTexParam("imgbuff", imgbuff->getTextureID());
		//s_texcontourProgram->SetGeometryTexParam("randMap", randTex);
		s_texcontourProgram->SetGeometryTexParam("ndoteMap", ndotetex->getTextureID());
		s_texcontourProgram->Enable();
		glBegin(GL_POINTS);
		for (int y = 0; y < vpH; y+=dsampleY)
		{
			int yt = (y + m_sclstep) % vpH;
			glVertex2f(0.f, (float)yt + 0.5);
		}
		glEnd();
		s_texcontourProgram->Disable();
		glDisable(GL_ALPHA_TEST);
	}

	//texturing contours vertically
	{
		glEnable(GL_ALPHA_TEST);
		glEnable(GL_BLEND);
		glBlendFunc(GL_DST_COLOR, GL_ZERO);
		glAlphaFunc(GL_LESS, 0.5f);
		glViewport(0, 0, vpW, vpH);
		s_texcontourVProgram->SetGeometryTexParam("positionMap", footprint->getTextureID());
		s_texcontourVProgram->SetGeometryTexParam("edgeMap", edgetex->getTextureID());
		s_texcontourVProgram->SetFragmentTexParam("strokeMap", contourstrokeTexture);
		s_texcontourVProgram->SetFragmentTexParam("depthMap", depthtex->getTextureID());
		s_texcontourVProgram->SetGeometryParam2f("scrSize", vpW, vpH);
		s_texcontourVProgram->SetGeometryParam1f("rsize", m_contourtexsize);
		s_texcontourVProgram->SetGeometryTexParam("imgbuff", imgbuff->getTextureID());
		//s_texcontourVProgram->SetGeometryTexParam("randMap", randTex);
		s_texcontourVProgram->SetGeometryTexParam("ndoteMap", ndotetex->getTextureID());
		s_texcontourVProgram->Enable();
		glBegin(GL_POINTS);
		for (int x = 0; x < vpW; x+=dsampleX)
		{
			int xt = (x + m_sclstep) % vpW;
			glVertex2f(0.f, (float)xt + 0.5);
		}
		glEnd();
		s_texcontourVProgram->Disable();
		glDisable(GL_ALPHA_TEST);
		glDisable(GL_BLEND);
	}

	m_sclstep++;
	
}

void FluidGridSimulator::NPRcoloring()
{
	//render contours
	//ContourRendering();
	glClearColor(1, 1, 1, 0.f); glClear(GL_COLOR_BUFFER_BIT);

	//////////////////////////////////////////////////////////
	//NPR colorize
	glEnable(GL_BLEND);
	glBlendFunc(GL_DST_COLOR, GL_ZERO);
	s_RenderQuadVertexProgram->Enable();
	s_toonshadeProgram->SetFragmentTexParam("imgbuff", imgbuff->getTextureID());
	s_toonshadeProgram->SetFragmentTexParam("depthMap", depthtex->getTextureID());
	s_toonshadeProgram->SetFragmentTexParam("colorband", nprcolorBand);
	s_toonshadeProgram->SetFragmentTexParam("ndoteMap", ndotetex->getTextureID());
	s_toonshadeProgram->SetFragmentParam3f("lightpos", lightPos[0], lightPos[1], lightPos[2]);
	s_toonshadeProgram->Enable();
	glViewport(0, 0, vpW, vpH);
	glBegin(GL_QUADS);
	glTexCoord2f(0,0); glVertex3f(-1,-1,0);
	glTexCoord2f(1,0); glVertex3f(1,-1,0);
	glTexCoord2f(1,1); glVertex3f(1,1,0);
	glTexCoord2f(0,1); glVertex3f(-1,1,0);		
	glEnd();
	s_toonshadeProgram->Disable();
	s_RenderQuadVertexProgram->Disable();
	glDisable(GL_BLEND);


	///////////////////////////////////////////////////////////////
	////make an initial point distribution and mark
	//{
	//	initpointtex->startRender();
	//	glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
	//	glViewport(0, 0, vpW, vpH);
	//	s_intertexturesProgram->SetVertexTexParam("positionMap", footprint->getTextureID());
	//	s_intertexturesProgram->SetVertexTexParam("imgbuff", imgbuff->getTextureID());
	//	s_intertexturesProgram->SetVertexParam2f("scrSize", vpW, vpH);
	//	s_intertexturesProgram->Enable();

	//	glCallList(samplingListId);

	//	s_intertexturesProgram->Disable();
	//	initpointtex->endRender();
	//}

	////texturing contours horizontally
	//{
	//	glEnable(GL_ALPHA_TEST);
	//	glAlphaFunc(GL_LESS, 0.1f);
	//	//glEnable(GL_BLEND);
	//	//glBlendFunc(GL_DST_COLOR, GL_ZERO);
	//	glViewport(0, 0, vpW, vpH);
	//	s_NPRcolorizeProgram->SetVertexTexParam("initpointMap", initpointtex->getTextureID());
	//	s_NPRcolorizeProgram->SetVertexTexParam("lightMap", imgbuff->getTextureID());
	//	s_NPRcolorizeProgram->SetVertexParam2f("scrSize", vpW, vpH);
	//	s_NPRcolorizeProgram->SetFragmentTexParam("splatMap", opaquediskTexture);
	//	s_NPRcolorizeProgram->SetFragmentTexParam("positionMap", footprint->getTextureID());
	//	s_NPRcolorizeProgram->Enable();
	//	glCallList(samplingListId);
	//	s_NPRcolorizeProgram->Disable();
	//	glDisable(GL_ALPHA_TEST);
	//	//glDisable(GL_BLEND);
	//}
}

void FluidGridSimulator::Stippling()
{
	//make an initial point distribution and mark
	{
		initpointtex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_intertexturesProgram->SetVertexTexParam("positionMap", footprint->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("imgbuff", imgbuff->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("ndoteMap", ndotetex->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("edgeMap", edgetex->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("toneMap", halftoneBand);
		s_intertexturesProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_intertexturesProgram->SetVertexParam1f("intensitymode", intensitymode);
		s_intertexturesProgram->SetVertexParam1f("a", gausssig);
		s_intertexturesProgram->SetVertexParam3f("lightpos", lightPos[0], lightPos[1], lightPos[2]);
		s_intertexturesProgram->SetVertexParam1f("CONTOUR_DENS_THRES", m_contourthres);
		s_intertexturesProgram->Enable();

		//glCallList(dispListId);					//need to be improved (i+=4; j+=4)
		glCallList(samplingListId);

		s_intertexturesProgram->Disable();
		initpointtex->endRender();
	}

	//mark covering region
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		marktex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glClearDepth(1.f);glClear(GL_DEPTH_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_markRegionProgram->SetVertexTexParam("initpointMap", initpointtex->getTextureID());
		s_markRegionProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_markRegionProgram->Enable();

		glCallList(samplingListId);					

		s_markRegionProgram->Disable();
		marktex->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}

	//advect 2nd view
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		tmpadvecttex_2v->startRender();
		glViewport(0, 0, vpW, vpH);
		glClearColor(0.f, 0.f, 0.f, 0.f); glClearDepth(1.f);
		glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
		s_advectSampPoint->SetGeometryTexParam("pointMap", pointCollection_2v->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("initpointMap", initpointtex->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("positionMap", footprint->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("velocityMap", velocity->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("offsetTexture", layerOffsetTexture);
		s_advectSampPoint->SetGeometryTexParam("markMap", marktex->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("densityMap", renderedDensity->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("imgbuff", imgbuff->getTextureID());
		s_advectSampPoint->SetGeometryParam3f("gridscale", 1.0 / m_size[0], 1.0 / m_size[1], 1.0 / m_size[2]);
		s_advectSampPoint->SetGeometryParam3f("gridcenter", m_pos[0], m_pos[1], m_pos[2]);
		s_advectSampPoint->SetGeometryParam3f("resolution", res_m, res_n, res_l);
		s_advectSampPoint->SetGeometryParam2f("scrSize", vpW, vpH);
		s_advectSampPoint->SetGeometryParam1f("d0", densityHaltonN);
		s_advectSampPoint->Enable();
		glCallList(dispListId);
		s_advectSampPoint->Disable();
		tmpadvecttex_2v->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}
	////main advect point collection
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		tmpadvecttex->startRender();
		glViewport(0, 0, vpW, vpH);
		glClearColor(0.f, 0.f, 0.f, 0.f); glClearDepth(1.f);
		glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
		s_advectSampPoint->SetGeometryTexParam("pointMap", pointCollection->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("initpointMap", initpointtex->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("positionMap", footprint->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("velocityMap", velocity->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("offsetTexture", layerOffsetTexture);
		s_advectSampPoint->SetGeometryTexParam("markMap", marktex->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("densityMap", renderedDensity->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("imgbuff", imgbuff->getTextureID());
		s_advectSampPoint->SetGeometryParam3f("gridscale", 1.0 / m_size[0], 1.0 / m_size[1], 1.0 / m_size[2]);
		s_advectSampPoint->SetGeometryParam3f("gridcenter", m_pos[0], m_pos[1], m_pos[2]);
		s_advectSampPoint->SetGeometryParam3f("resolution", res_m, res_n, res_l);
		s_advectSampPoint->SetGeometryParam2f("scrSize", vpW, vpH);
		s_advectSampPoint->SetGeometryParam1f("d0", densityHaltonN);
		s_advectSampPoint->Enable();
		glCallList(dispListId);
		s_advectSampPoint->Disable();
		tmpadvecttex->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}

	//distribute points 2nd
	{
		advectCollection_2v->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		s_distPointProgram2v->SetVertexTexParam("initpointMap", tmpadvecttex_2v->getTextureID());
		s_distPointProgram2v->SetVertexParam2f("scrSize", vpW, vpH);
		s_distPointProgram2v->SetVertexTexParam("markMap", marktex->getTextureID());
		s_distPointProgram2v->SetVertexTexParam("imgbuff", imgbuff->getTextureID());
		s_distPointProgram2v->SetFragmentParam1f("para", 1.f);
		s_distPointProgram2v->SetFragmentParam1f("d0", densityHaltonN);
		s_distPointProgram2v->Enable();
		glCallList(dispListId);
		s_distPointProgram2v->Disable();
		advectCollection_2v->endRender();
	}

	////main distribute points
	{
		advectCollection->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		s_distPointProgram->SetVertexTexParam("pointMap", tmpadvecttex->getTextureID());
		s_distPointProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_distPointProgram->SetVertexTexParam("markMap", marktex->getTextureID());
		s_distPointProgram->SetVertexTexParam("imgbuff", imgbuff->getTextureID());
		s_distPointProgram->SetFragmentParam1f("d0", densityHaltonN);
		s_distPointProgram->Enable();
		glCallList(dispListId);
		s_distPointProgram->Disable();
		advectCollection->endRender();
	}


	//////////////////////////////////////////////////////////////////////////
	//main second viewport
	glClearColor(1.f, 1.f, 1.f, 1.f); glClear(GL_COLOR_BUFFER_BIT);
	glEnable(GL_POINT_SMOOTH);
	glPointSize(1.f);
	glEnable(GL_BLEND);
	glBlendFunc(GL_DST_COLOR, GL_ZERO);
	glViewport(0, 0, vpW, vpH);
	s_stippleProgram->SetVertexTexParam("pointMap", advectCollection_2v->getTextureID());
	s_stippleProgram->SetVertexParam2f("scrSize", vpW, vpH);
	s_stippleProgram->Enable();
	glCallList(dispListId);
	s_stippleProgram->Disable();
	glDisable(GL_BLEND);
	glDisable(GL_POINT_SMOOTH);
	glPointSize(1.f);

	//the main viewport
	glEnable(GL_POINT_SMOOTH);
	glPointSize(2.f);
	glEnable(GL_BLEND);
	glBlendFunc(GL_DST_COLOR, GL_ZERO);
	glViewport(vpW, 0, vpW, vpH);
	s_stippleProgram->SetVertexTexParam("pointMap", advectCollection->getTextureID());
	s_stippleProgram->SetVertexParam2f("scrSize", vpW, vpH);
	s_stippleProgram->Enable();
	glCallList(dispListId);
	s_stippleProgram->Disable();
	glDisable(GL_BLEND);
	glDisable(GL_POINT_SMOOTH);
	glPointSize(1.f);

	//////////////////////////////////////////////////////////////////////////

	//renderRealSmoke(vpW, vpH, NULL, vpW, 0);

	//swap
	RenderTexture* tmp = advectCollection;
	advectCollection = pointCollection;
	pointCollection = tmp;

	//swap 2nd
	tmp = advectCollection_2v;
	advectCollection_2v = pointCollection_2v;
	pointCollection_2v = tmp;
}

void FluidGridSimulator::Fading()
{
}

void FluidGridSimulator::SolidContour()
{
	//make an initial point distribution and mark
	{
		initpointtex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_intertexturesProgram->SetVertexTexParam("positionMap", footprint->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("imgbuff", imgbuff->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("ndoteMap", ndotetex->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("edgeMap", edgetex->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("toneMap", halftoneBand);
		s_intertexturesProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_intertexturesProgram->SetVertexParam1f("intensitymode", intensitymode);
		s_intertexturesProgram->SetVertexParam1f("a", gausssig);
		s_intertexturesProgram->Enable();

		//glCallList(dispListId);					//need to be improved (i+=4; j+=4)
		glCallList(samplingListId);

		s_intertexturesProgram->Disable();
		initpointtex->endRender();
	}

	//mark covering region
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		marktex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glClearDepth(1.f);glClear(GL_DEPTH_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_markRegionProgram->SetVertexTexParam("initpointMap", initpointtex->getTextureID());
		s_markRegionProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_markRegionProgram->Enable();

		//glCallList(dispListId);
		glCallList(samplingListId);					

		s_markRegionProgram->Disable();
		marktex->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}

	//advect point collection
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		tmpadvecttex->startRender();
		glViewport(0, 0, vpW, vpH);
		glClearColor(0.f, 0.f, 0.f, 0.f); glClearDepth(1.f);
		glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
		s_advectSampPoint->SetGeometryTexParam("pointMap", pointCollection->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("initpointMap", initpointtex->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("positionMap", footprint->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("velocityMap", velocity->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("offsetTexture", layerOffsetTexture);
		s_advectSampPoint->SetGeometryTexParam("markMap", marktex->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("densityMap", renderedDensity->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("imgbuff", imgbuff->getTextureID());
		s_advectSampPoint->SetGeometryParam3f("gridscale", 1.0 / m_size[0], 1.0 / m_size[1], 1.0 / m_size[2]);
		s_advectSampPoint->SetGeometryParam3f("gridcenter", m_pos[0], m_pos[1], m_pos[2]);
		s_advectSampPoint->SetGeometryParam3f("resolution", res_m, res_n, res_l);
		s_advectSampPoint->SetGeometryParam2f("scrSize", vpW, vpH);
		s_advectSampPoint->SetGeometryParam1f("d0", densityHaltonN);
		s_advectSampPoint->Enable();
		glCallList(dispListId);
		s_advectSampPoint->Disable();
		tmpadvecttex->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}

	//distribute points
	{
		advectCollection->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		s_distPointProgram->SetVertexTexParam("pointMap", tmpadvecttex->getTextureID());
		s_distPointProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_distPointProgram->SetVertexTexParam("markMap", marktex->getTextureID());
		s_distPointProgram->SetVertexTexParam("imgbuff", imgbuff->getTextureID());
		s_distPointProgram->SetFragmentParam1f("d0", densityHaltonN);
		s_distPointProgram->Enable();
		glCallList(dispListId);
		s_distPointProgram->Disable();
		advectCollection->endRender();
	}

	//swap
	RenderTexture* tmp = advectCollection;
	advectCollection = pointCollection;
	pointCollection = tmp;

	//render with billboards
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glClearDepth(1.f);
		glClearColor(0.69f, 0.8f, 0.9f, 0.f);//(1.f, 1.f, 1.f, 0.f); 
		glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
		glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
		glEnable(GL_ALPHA_TEST);
		glAlphaFunc(GL_GREATER, 0.f);
		glViewport(0, 0, vpW, vpH);
		s_solidcontourProgram->SetVertexTexParam("positionMap", footprint->getTextureID());
		s_solidcontourProgram->SetVertexTexParam("pointMap", pointCollection->getTextureID());
		s_solidcontourProgram->SetVertexParam1f("recsize", recsize);
		s_solidcontourProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_solidcontourProgram->SetVertexParam3f("gridscale", 1.0 / m_size[0], 1.0 / m_size[1], 1.0 / m_size[2]);
		s_solidcontourProgram->SetVertexParam3f("gridcenter", m_pos[0], m_pos[1], m_pos[2]);
		s_solidcontourProgram->SetVertexParam3f("resolution", res_m, res_n, res_l);
		s_solidcontourProgram->SetVertexTexParam("densityMap", renderedDensity->getTextureID());
		s_solidcontourProgram->SetVertexTexParam("offsetTexture", layerOffsetTexture);
		s_solidcontourProgram->SetVertexTexParam("randMap", randTex);
		s_solidcontourProgram->SetFragmentTexParam("billboard", billBoard->getTextureID());
		s_solidcontourProgram->SetFragmentParam3f("lightpos", lightPos[0], lightPos[1], lightPos[2]);
		s_solidcontourProgram->Enable();

		glCallList(dispListId);

		s_solidcontourProgram->Disable();
		glDisable(GL_ALPHA_TEST);
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}
}

void FluidGridSimulator::Hatching()
{
/*
	//edge detection filtering
	{
		edgetex->startRender();
		s_RenderQuadVertexProgram->Enable();
		s_edgedetectProgram->SetFragmentTexParam("ndotetext", ndotetex->getTextureID());
		s_edgedetectProgram->SetFragmentTexParam("imgbuff", imgbuff->getTextureID());
		s_edgedetectProgram->SetFragmentParam2f("scrSize", vpW, vpH);
		s_edgedetectProgram->SetFragmentParam1f("edgemode", 0);
		s_edgedetectProgram->Enable();
		glViewport(0, 0, vpW, vpH);
		glBegin(GL_QUADS);
		glTexCoord2f(0,0); glVertex3f(-1,-1,0);
		glTexCoord2f(1,0); glVertex3f(1,-1,0);
		glTexCoord2f(1,1); glVertex3f(1,1,0);
		glTexCoord2f(0,1); glVertex3f(-1,1,0);		
		glEnd();
		s_edgedetectProgram->Disable();
		s_RenderQuadVertexProgram->Disable();
		edgetex->endRender();
	}

	if (edgemode == 1)
		//combine edges
	{
		edgetex->startRender();
		s_RenderQuadVertexProgram->Enable();
		s_edgedetectProgram->SetFragmentTexParam("ndotetext", ndotetex->getTextureID());
		s_edgedetectProgram->SetFragmentTexParam("imgbuff", imgbuff->getTextureID());
		s_edgedetectProgram->SetFragmentParam2f("scrSize", vpW, vpH);
		s_edgedetectProgram->SetFragmentParam1f("edgemode", 1);
		s_edgedetectProgram->SetFragmentParam1f("joinmode", joinmode);
		s_edgedetectProgram->SetFragmentTexParam("edgeMap", edgetex->getTextureID());
		s_edgedetectProgram->SetFragmentParam1f("a", gausssig);
		s_edgedetectProgram->Enable();
		glViewport(0, 0, vpW, vpH);
		glBegin(GL_QUADS);
		glTexCoord2f(0,0); glVertex3f(-1,-1,0);
		glTexCoord2f(1,0); glVertex3f(1,-1,0);
		glTexCoord2f(1,1); glVertex3f(1,1,0);
		glTexCoord2f(0,1); glVertex3f(-1,1,0);		
		glEnd();
		s_edgedetectProgram->Disable();
		s_RenderQuadVertexProgram->Disable();
		edgetex->endRender();
	}
*/
	ContourRendering();
	////////////////////hatching
	//make an initial point distribution and mark
	{
		initpointtex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_intertexturesProgram->SetVertexTexParam("positionMap", footprint->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("imgbuff", imgbuff->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("ndoteMap", ndotetex->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("edgeMap", edgetex->getTextureID());
		s_intertexturesProgram->SetVertexTexParam("toneMap", halftoneBand);
		s_intertexturesProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_intertexturesProgram->SetVertexParam1f("intensitymode", intensitymode);
		s_intertexturesProgram->SetVertexParam1f("a", gausssig);
		s_intertexturesProgram->SetVertexParam1f("CONTOUR_DENS_THRES", m_contourthres);
		s_intertexturesProgram->Enable();

		//glCallList(dispListId);					//need to be improved (i+=4; j+=4)
		glCallList(samplingListId);

		s_intertexturesProgram->Disable();
		initpointtex->endRender();
	}

	//mark covering region
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		marktex->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		glClearDepth(1.f);glClear(GL_DEPTH_BUFFER_BIT);
		glViewport(0, 0, vpW, vpH);
		s_markRegionProgram->SetVertexTexParam("initpointMap", initpointtex->getTextureID());
		s_markRegionProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_markRegionProgram->Enable();

		//glCallList(dispListId);
		glCallList(samplingListId);					

		s_markRegionProgram->Disable();
		marktex->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}

	//advect point collection
	{
		glEnable(GL_DEPTH_TEST); glDepthMask(GL_TRUE);
		glDisable(GL_BLEND);
		tmpadvecttex->startRender();
		glViewport(0, 0, vpW, vpH);
		glClearColor(0.f, 0.f, 0.f, 0.f); glClearDepth(1.f);
		glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
		s_advectSampPoint->SetGeometryTexParam("pointMap", pointCollection->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("initpointMap", initpointtex->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("positionMap", footprint->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("velocityMap", velocity->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("offsetTexture", layerOffsetTexture);
		s_advectSampPoint->SetGeometryTexParam("markMap", marktex->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("densityMap", renderedDensity->getTextureID());
		s_advectSampPoint->SetGeometryTexParam("imgbuff", imgbuff->getTextureID());
		s_advectSampPoint->SetGeometryParam3f("gridscale", 1.0 / m_size[0], 1.0 / m_size[1], 1.0 / m_size[2]);
		s_advectSampPoint->SetGeometryParam3f("gridcenter", m_pos[0], m_pos[1], m_pos[2]);
		s_advectSampPoint->SetGeometryParam3f("resolution", res_m, res_n, res_l);
		s_advectSampPoint->SetGeometryParam2f("scrSize", vpW, vpH);
		s_advectSampPoint->SetGeometryParam1f("d0", densityHaltonN);
		s_advectSampPoint->Enable();
		glCallList(dispListId);
		s_advectSampPoint->Disable();
		tmpadvecttex->endRender();
		glDisable(GL_DEPTH_TEST); glDepthMask(GL_FALSE);
	}

	//distribute points
	{
		advectCollection->startRender();
		glClearColor(0.f, 0.f, 0.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		s_distPointProgram->SetVertexTexParam("pointMap", tmpadvecttex->getTextureID());
		s_distPointProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_distPointProgram->SetVertexTexParam("markMap", marktex->getTextureID());
		s_distPointProgram->SetVertexTexParam("imgbuff", imgbuff->getTextureID());
		s_distPointProgram->SetFragmentParam1f("d0", densityHaltonN);
		s_distPointProgram->Enable();
		glCallList(dispListId);
		s_distPointProgram->Disable();
		advectCollection->endRender();
	}

	//swap
	RenderTexture* tmp = advectCollection;
	advectCollection = pointCollection;
	pointCollection = tmp;

	//hatching with by accessing the texture stored curvature info
	{
		//glClearColor(0.2f, 0.2f, 0.2f, 0.2f); glClear(GL_COLOR_BUFFER_BIT);
		//glClearColor(0.95f, 0.77f, 0.52f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		//glClearColor(1.f, 1.f, 1.f, 0.f); glClear(GL_COLOR_BUFFER_BIT);
		//glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
		glEnable(GL_BLEND);
		glBlendFunc(GL_DST_COLOR, GL_ZERO);
		//glBlendFunc(GL_DST_COLOR, GL_ONE_MINUS_DST_COLOR);
		glEnable(GL_ALPHA_TEST);
		glAlphaFunc(GL_LESS, 0.3f);
		glViewport(0, 0, vpW, vpH);
		s_hatchingProgram->SetVertexTexParam("positionMap", footprint->getTextureID());
		s_hatchingProgram->SetVertexTexParam("pointMap", pointCollection->getTextureID());
		s_hatchingProgram->SetVertexTexParam("densityMap", renderedDensity->getTextureID());
		s_hatchingProgram->SetVertexTexParam("ndoteMap", ndotetex->getTextureID());
		s_hatchingProgram->SetVertexTexParam("offsetTexture", layerOffsetTexture);
		s_hatchingProgram->SetVertexTexParam("hesrow1", hessian_r1->getTextureID());
		s_hatchingProgram->SetVertexTexParam("hesrow2", hessian_r2->getTextureID());
		s_hatchingProgram->SetVertexTexParam("hesrow3", hessian_r3->getTextureID());
		//s_hatchingProgram->SetVertexTexParam("imgbuff", imgbuff->getTextureID());
		s_hatchingProgram->SetVertexParam3f("gridscale", 1.0 / m_size[0], 1.0 / m_size[1], 1.0 / m_size[2]);
		s_hatchingProgram->SetVertexParam3f("gridcenter", m_pos[0], m_pos[1], m_pos[2]);
		s_hatchingProgram->SetVertexParam3f("resolution", res_m, res_n, res_l);
		s_hatchingProgram->SetVertexParam1f("recsize", recsize);
		s_hatchingProgram->SetVertexParam1f("crsize", m_contourtexsize);
		s_hatchingProgram->SetVertexParam2f("scrSize", vpW, vpH);
		s_hatchingProgram->SetVertexParam1f("hatchmode", hatchmode);
		//s_hatchingProgram->SetVertexParam1f("curvaturemode", detectMode);
		s_hatchingProgram->SetVertexTexParam("edgeMap", edgetex->getTextureID());
		s_hatchingProgram->SetFragmentTexParam("strokeMap", strokeTexture);
		//s_hatchingProgram->SetGeometryTexParam("randMap", randTex);
		//s_hatchingProgram->SetFragmentTexParam("colorMap", colorstrokeBand);
		//s_hatchingProgram->SetFragmentTexParam("depthMap", depthtex->getTextureID());
		//s_hatchingProgram->SetFragmentTexParam("lightMap", imgbuff->getTextureID());
		//s_hatchingProgram->SetFragmentTexParam("ndoteMap", ndotetex->getTextureID());
		//s_hatchingProgram->SetFragmentParam2f("scrSize", vpW, vpH);
		s_hatchingProgram->SetVertexParam1f("CONTOUR_DENS_THRES", m_contourthres);
		s_hatchingProgram->SetVertexParam1f("cellsize", sizeX / res_m);
		s_hatchingProgram->SetVertexParam1f("stroketype", innerStrokeType);
		s_hatchingProgram->Enable();

		glCallList(dispListId);

		s_hatchingProgram->Disable();
		glDisable(GL_BLEND);
		glDisable(GL_ALPHA_TEST);
	}

/*
	//texturing contours horizontally
	{
		glEnable(GL_ALPHA_TEST);
		glAlphaFunc(GL_LESS, 0.5f);
		glViewport(0, 0, vpW, vpH);
		s_texcontourProgram->SetGeometryTexParam("positionMap", footprint->getTextureID());
		s_texcontourProgram->SetGeometryTexParam("edgeMap", edgetex->getTextureID());
		s_texcontourProgram->SetFragmentTexParam("strokeMap", contourstrokeTexture);
		s_texcontourProgram->SetFragmentTexParam("depthMap", depthtex->getTextureID());
		s_texcontourProgram->SetGeometryParam2f("scrSize", vpW, vpH);
		s_texcontourProgram->SetGeometryParam1f("rsize", dsampleY);
		s_texcontourProgram->SetGeometryTexParam("imgbuff", imgbuff->getTextureID());
		s_texcontourProgram->Enable();
		glBegin(GL_POINTS);
		for (int y = 0; y < vpH; y+=dsampleY)
			glVertex2f(0.f, (float)y + 0.5);
		glEnd();
		s_texcontourProgram->Disable();
		glDisable(GL_ALPHA_TEST);
	}

	//texturing contours vertically
	{
		glEnable(GL_ALPHA_TEST);
		glAlphaFunc(GL_LESS, 0.5f);
		glViewport(0, 0, vpW, vpH);
		s_texcontourVProgram->SetGeometryTexParam("positionMap", footprint->getTextureID());
		s_texcontourVProgram->SetGeometryTexParam("edgeMap", edgetex->getTextureID());
		s_texcontourVProgram->SetFragmentTexParam("strokeMap", contourstrokeTexture);
		s_texcontourVProgram->SetFragmentTexParam("depthMap", depthtex->getTextureID());
		s_texcontourVProgram->SetGeometryParam2f("scrSize", vpW, vpH);
		s_texcontourVProgram->SetGeometryParam1f("rsize", dsampleX);
		s_texcontourVProgram->SetGeometryTexParam("imgbuff", imgbuff->getTextureID());
		s_texcontourVProgram->Enable();
		glBegin(GL_POINTS);
		for (int x = 0; x < vpW; x+=dsampleX)
			glVertex2f(0.f, (float)x + 0.5);
		glEnd();
		s_texcontourVProgram->Disable();
		glDisable(GL_ALPHA_TEST);
	}
*/
}

void FluidGridSimulator::renderGridPoint(int w, int h)
{
	s_RenderQuadVertexProgram->Enable();
	//compute density gradient
	density->startRender();
	s_denGraProgram->SetFragmentTexParam("densityMap", density->getTextureID());
	s_denGraProgram->SetFragmentTexParam("colorbandMap", strokeTexture);
	s_denGraProgram->Enable();
	s_denGraProgram->EnableFragmentTexParam("densityMap");
	s_denGraProgram->EnableFragmentTexParam("colorbandMap");
	for(int l = 0; l < res_l; l++)
	{
		gridLayers[l]->DrawQuad();
	}
	s_denGraProgram->DisableFragmentTexParam("densityMap");
	s_denGraProgram->DisableFragmentTexParam("colorbandMap");
	s_denGraProgram->Disable();
	density->endRender();
	s_RenderQuadVertexProgram->Disable();

	//render bounding cube
	if(renderBoundingCube)
	{
		glViewport(0, 0, w, h);
		glClear(GL_COLOR_BUFFER_BIT);
		glColor4f(0,0.8,0.7,0.2);
		glPushMatrix();
		glTranslatef(posX, posZ, posY);
		glScalef(sizeX, sizeZ, sizeY);
		glutWireCube(1);
		glPopMatrix();
	}

	s_gridpointProgram->Enable();
	s_gridpointProgram->SetFragmentTexParam("density", density->getTextureID());
	s_gridpointProgram->SetFragmentTexParam("offsetTexture", layerOffsetTexture);
	s_gridpointProgram->SetFragmentTexParam("colorbandMap", strokeTexture);
	s_gridpointProgram->EnableFragmentTexParam("density");
	s_gridpointProgram->EnableFragmentTexParam("offsetTexture");
	s_gridpointProgram->EnableFragmentTexParam("colorbandMap");
	glViewport(0, 0, w, h);
	glEnable(GL_BLEND);
	glBlendFunc(GL_ONE, GL_ONE);
	glBegin(GL_POINTS);
	float d = sizeX / res_m;
	float x = -sizeX / 2.f + d / 2.f;
	float z = -sizeY / 2.f + d / 2.f;
	float y = -sizeZ / 2.f + d / 2.f;
	for (int i = 0; i < res_m; i++)
		for (int j = 0; j < res_l; j++)
			for (int k = 0; k < res_n ; k++)
			{
				glTexCoord3f(i + 0.5, k + 0.5, j + 1);
				glVertex3f(x + i * d, y + j * d, z + k * d);
			}
	glEnd();
	glDisable(GL_BLEND);
	s_gridpointProgram->DisableFragmentTexParam("density");
	s_gridpointProgram->DisableFragmentTexParam("offsetTexture");
	s_gridpointProgram->DisableFragmentTexParam("colorbandMap");
	s_gridpointProgram->Disable();
}

void FluidGridSimulator::resetTemperature()
{
	temperature->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);	
	temperature->endRender();

	temperature->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);	
	temperature->endRender();
}

void FluidGridSimulator::resetDensity()
{
	density->reset();
	initializeImpulse();
	m_time = 0.f;
}

void FluidGridSimulator::resetVelocity()
{
	velocity->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);	
	velocity->endRender();

	velocity->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);	
	velocity->endRender();

	pressure->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);
	pressure->endRender();

	pressure->startRender();
	glClearColor(0,0,0,0);
	glClear(GL_COLOR_BUFFER_BIT);
	pressure->endRender();
}

//the color band store color used to set color for grid point in gradient rendering
void FluidGridSimulator::createColorBand()
{
	char silhouetteTexFile[100];
	char shadeTexFile[100];
	switch (strokeTexType)
	{
	case 0:
		sprintf(silhouetteTexFile,"contourstroke_8.jpg");
		break;
	case 1:
		sprintf(silhouetteTexFile,"contourstroke_9.jpg");
		break;
	case 2:
		sprintf(silhouetteTexFile,"contourstroke_10_1.jpg");
		break;
	}
	sprintf(shadeTexFile,"stroke_3.bmp");

	ILuint texid;
	ILboolean success;
	ilInit();
	ilGenImages(1, &texid);
	ilBindImage(texid);
	success = ilLoadImage(shadeTexFile);
	if (success)
	{
		success = ilConvertImage(IL_RGB, IL_FLOAT);
		if (!success)
			printf("Can not convert stroke image");
		glGenTextures(1, &strokeTexture);
		glBindTexture(GL_TEXTURE_2D, strokeTexture);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
		glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, ilGetInteger(IL_IMAGE_WIDTH),
			ilGetInteger(IL_IMAGE_HEIGHT), 0, ilGetInteger(IL_IMAGE_FORMAT), GL_FLOAT,
			ilGetData());
		glBindTexture(GL_TEXTURE_2D, 0);
	}
	else
	{		
		printf("Can not load stroke image\n");
	}
	ilDeleteImages(1, &texid);

	//texture of stroke for contours
	ilInit();
	ilGenImages(1, &texid);
	ilBindImage(texid);
	success = ilLoadImage(silhouetteTexFile);
	if (success)
	{
		success = ilConvertImage(IL_RGB, IL_FLOAT);
		if (!success)
			printf("Can not convert contour stroke image");
		glGenTextures(1, &contourstrokeTexture);
		glBindTexture(GL_TEXTURE_2D, contourstrokeTexture);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
		glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, ilGetInteger(IL_IMAGE_WIDTH),
			ilGetInteger(IL_IMAGE_HEIGHT), 0, ilGetInteger(IL_IMAGE_FORMAT), GL_FLOAT,
			ilGetData());
		glBindTexture(GL_TEXTURE_2D, 0);
	}
	else
	{		
		printf("Can not load contour stroke image\n");
	}
	ilDeleteImages(1, &texid);

	//texture for splats
	ilInit();
	ilGenImages(1, &texid);
	ilBindImage(texid);
	success = ilLoadImage("splat.bmp");
	if (success)
	{
		success = ilConvertImage(IL_RGB, IL_FLOAT);
		if (!success)
			printf("Can not convert color band image");
		glGenTextures(1, &splatTexture);
		glBindTexture(GL_TEXTURE_2D, splatTexture);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
		glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, ilGetInteger(IL_IMAGE_WIDTH),
			ilGetInteger(IL_IMAGE_HEIGHT), 0, ilGetInteger(IL_IMAGE_FORMAT), GL_FLOAT,
			ilGetData());
		glBindTexture(GL_TEXTURE_2D, 0);
	}
	else
	{		
		printf("Can not load splat image\n");
	}
	ilDeleteImages(1, &texid);

	//texture for opaque disks
	ilInit();
	ilGenImages(1, &texid);
	ilBindImage(texid);
	success = ilLoadImage("opaquedisk.bmp");
	if (success)
	{
		success = ilConvertImage(IL_RGB, IL_FLOAT);
		if (!success)
			printf("Can not convert color band image");
		glGenTextures(1, &opaquediskTexture);
		glBindTexture(GL_TEXTURE_2D, opaquediskTexture);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
		glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, ilGetInteger(IL_IMAGE_WIDTH),
			ilGetInteger(IL_IMAGE_HEIGHT), 0, ilGetInteger(IL_IMAGE_FORMAT), GL_FLOAT,
			ilGetData());
		glBindTexture(GL_TEXTURE_2D, 0);
	}
	else
	{		
		printf("Can not load opaque disk image\n");
	}
	ilDeleteImages(1, &texid);

	ilInit();
	ilGenImages(1, &texid);
	ilBindImage(texid);
	success = ilLoadImage("nprband_small.bmp");
	if (success)
	{
		success = ilConvertImage(IL_RGB, IL_FLOAT);
		if (!success)
			printf("Can not convert color band image");
		glGenTextures(1, &nprcolorBand);
		glBindTexture(GL_TEXTURE_2D, nprcolorBand);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
		glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, ilGetInteger(IL_IMAGE_WIDTH),
			ilGetInteger(IL_IMAGE_HEIGHT), 0, ilGetInteger(IL_IMAGE_FORMAT), GL_FLOAT,
			ilGetData());
		glBindTexture(GL_TEXTURE_2D, 0);
	}
	else
	{		
		printf("Can not load npr band image\n");
	}
	ilDeleteImages(1, &texid);

	//color band for transfer function
	ilInit();
	ilGenImages(1, &texid);
	ilBindImage(texid);
	success = ilLoadImage("halftone.jpg");
	if (success)
	{
		success = ilConvertImage(IL_RGB, IL_FLOAT);
		if (!success)
			printf("Can not convert half tone band image");
		glGenTextures(1, &halftoneBand);
		glBindTexture(GL_TEXTURE_2D, halftoneBand);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
		glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, ilGetInteger(IL_IMAGE_WIDTH),
			ilGetInteger(IL_IMAGE_HEIGHT), 0, ilGetInteger(IL_IMAGE_FORMAT), GL_FLOAT,
			ilGetData());
		glBindTexture(GL_TEXTURE_2D, 0);
	}
	else
	{		
		printf("Can not load half tone image\n");
	}
	ilDeleteImages(1, &texid);

	//color band for coloring strokes
	ilInit();
	ilGenImages(1, &texid);
	ilBindImage(texid);
	success = ilLoadImage("strokecolor.jpg");
	if (success)
	{
		success = ilConvertImage(IL_RGB, IL_FLOAT);
		if (!success)
			printf("Can not convert half tone band image");
		glGenTextures(1, &colorstrokeBand);
		glBindTexture(GL_TEXTURE_2D, colorstrokeBand);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
		glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
		glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, ilGetInteger(IL_IMAGE_WIDTH),
			ilGetInteger(IL_IMAGE_HEIGHT), 0, ilGetInteger(IL_IMAGE_FORMAT), GL_FLOAT,
			ilGetData());
		glBindTexture(GL_TEXTURE_2D, 0);
	}
	else
	{		
		printf("Can not load half tone image\n");
	}
	ilDeleteImages(1, &texid);
}

// the offset layer texture stores the texture coordinates in the flat3D texture of each grid layer origo
// this is a helper texture to easily find an X,Y coordinate in the flat3D texture of a point (x,y,z) in the 3D grid:
// (X,Y) = tex(layerOffsetTexture, z).xy + (x,y)
void FluidGridSimulator::createLayerOffsetTexture()
{
	glGenTextures(1, &layerOffsetTexture);

	glBindTexture(GL_TEXTURE_RECTANGLE_EXT, layerOffsetTexture);

	int M = res_m + 2; //2 pixel bound edges
	int N = res_n + 2; //2 pixel bound edges
	int L = res_l + 2; //2 bound layers

	float* data = new float[L * 4];

	int x = 0;
	int y = 0;
	for(int i = 0; i < L; i++)
	{
		data[i * 4] = x * M + 1;
		data[i * 4 + 1] = y * N + 1;
		data[i * 4 + 2] = 0;
		data[i * 4 + 3] = 0;

		x++;
		if(x == numTilesX)
		{
			x = 0;
			y++;
		}
	}

	glTexImage2D(GL_TEXTURE_RECTANGLE_EXT, 0, GL_RGBA16F_ARB, L, 1, 0, GL_RGBA, GL_FLOAT, data);
	glTexParameterf(GL_TEXTURE_RECTANGLE_EXT, GL_TEXTURE_WRAP_S, GL_CLAMP);
	glTexParameterf(GL_TEXTURE_RECTANGLE_EXT, GL_TEXTURE_WRAP_T, GL_CLAMP);
	glTexParameterf(GL_TEXTURE_RECTANGLE_EXT, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
	glTexParameterf(GL_TEXTURE_RECTANGLE_EXT, GL_TEXTURE_MIN_FILTER, GL_NEAREST);

	glBindTexture(GL_TEXTURE_RECTANGLE_EXT, 0);
}

//return the velocity, density and pressure value of the upmost inner layer
//used in distributed simulation. This layer will be set as a lower boundary layer to an adjacent fluid simulator
void FluidGridSimulator::getUpperLayer(void* data)
{
	char* pixeldata = (char*) data;
	int x = gridLayers[res_l - 1]->columnNum * (res_m + 2) + 1;
	int y = gridLayers[res_l - 1]->rowNum * (res_n + 2) + 1;

	velocity->startCopy();
	glReadPixels(x, y, res_m, res_n, GL_RGB, GL_HALF_FLOAT_ARB, pixeldata);
	GLenum error = glGetError();
	velocity->endCopy();

	density->startCopy();
	glReadPixels(x, y, res_m, res_n, GL_RED, GL_HALF_FLOAT_ARB, &pixeldata[res_m*res_n*3*2]);
	error = glGetError();
	density->endCopy();	

	pressure->startCopy();
	glReadPixels(x, y, res_m, res_n, GL_RED, GL_HALF_FLOAT_ARB, &pixeldata[res_m*res_n*4*2]);
	error = glGetError();
	pressure->endCopy();	
}

//return the velocity, density and pressure value of the lowest inner layer
//used in distributed simulation. This layer will be set as an upper boundary layer to an adjacent fluid simulator
void FluidGridSimulator::getLowerLayer(void* data)
{
	char* pixeldata = (char*) data;
	int x = res_m + 3;
	int y = 1;

	velocity->startCopy();
	glReadPixels(x, y, res_m, res_n, GL_RGB, GL_HALF_FLOAT_ARB, pixeldata);
	GLenum error = glGetError();
	velocity->endCopy();

	density->startCopy();
	glReadPixels(x, y, res_m, res_n, GL_RED, GL_HALF_FLOAT_ARB, &pixeldata[res_m*res_n*3*2]);
	error = glGetError();
	density->endCopy();

	pressure->startCopy();
	glReadPixels(x, y, res_m, res_n, GL_RED, GL_HALF_FLOAT_ARB, &pixeldata[res_m*res_n*4*2]);
	error = glGetError();
	pressure->endCopy();
}

//set the velocity, density and pressure value of the upper boundary layer
//used in distributed simulation. see getLowerLayer()
void FluidGridSimulator::setUpperBoundLayer(void* pixeldata)
{
	glPixelTransferf(GL_ALPHA_SCALE, 0);
	int x = boundLayers[0]->columnNum * (res_m + 2) + 1;
	int y = boundLayers[0]->rowNum * (res_n + 2) + 1;

	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, velocity->getRenderTextureID());
	glTexSubImage2D(GL_TEXTURE_RECTANGLE_ARB, 0, x, y, res_m, res_n, GL_RGB, GL_HALF_FLOAT_ARB, pixeldata);
	GLenum error = glGetError();
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, 0);

	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, velocity->getTextureID());
	glTexSubImage2D(GL_TEXTURE_RECTANGLE_ARB, 0, x, y, res_m, res_n, GL_RGB, GL_HALF_FLOAT_ARB, pixeldata);
	error = glGetError();
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, 0);

	char* data = (char*) pixeldata;
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, density->getTextureID());
	glTexSubImage2D(GL_TEXTURE_RECTANGLE_ARB, 0, x, y, res_m, res_n, GL_RED, GL_HALF_FLOAT_ARB, &data[res_m*res_n*3*2]);
	error = glGetError();
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, 0);

	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, density->getRenderTextureID());
	glTexSubImage2D(GL_TEXTURE_RECTANGLE_ARB, 0, x, y, res_m, res_n, GL_RED, GL_HALF_FLOAT_ARB, &data[res_m*res_n*3*2]);
	error = glGetError();
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, 0);

	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, pressure->getTextureID());
	glTexSubImage2D(GL_TEXTURE_RECTANGLE_ARB, 0, x, y, res_m, res_n, GL_RED, GL_HALF_FLOAT_ARB, &data[res_m*res_n*4*2]);
	error = glGetError();
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, 0);

	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, pressure->getRenderTextureID());
	glTexSubImage2D(GL_TEXTURE_RECTANGLE_ARB, 0, x, y, res_m, res_n, GL_RED, GL_HALF_FLOAT_ARB, &data[res_m*res_n*4*2]);
	error = glGetError();
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, 0);
	glPixelTransferf(GL_ALPHA_SCALE, 1);
}

//set the velocity, density and pressure value of the lower boundary layer
//used in distributed simulation. see getUpperLayer()
void FluidGridSimulator::setLowerBoundLayer(void* pixeldata)
{
	glPixelTransferf(GL_ALPHA_SCALE, 0);
	int x = 1;
	int y = 1;

	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, velocity->getRenderTextureID());
	glTexSubImage2D(GL_TEXTURE_RECTANGLE_ARB, 0, x, y, res_m, res_n, GL_RGB, GL_HALF_FLOAT_ARB, pixeldata);
	GLenum error = glGetError();
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, 0);

	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, velocity->getTextureID());
	glTexSubImage2D(GL_TEXTURE_RECTANGLE_ARB, 0, x, y, res_m, res_n, GL_RGB, GL_HALF_FLOAT_ARB, pixeldata);
	error = glGetError();
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, 0);

	char* data = (char*) pixeldata;

	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, density->getTextureID());
	glTexSubImage2D(GL_TEXTURE_RECTANGLE_ARB, 0, x, y, res_m, res_n, GL_RED, GL_HALF_FLOAT_ARB, &data[res_m*res_n*3*2]);
	error = glGetError();
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, 0);

	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, density->getRenderTextureID());
	glTexSubImage2D(GL_TEXTURE_RECTANGLE_ARB, 0, x, y, res_m, res_n, GL_RED, GL_HALF_FLOAT_ARB, &data[res_m*res_n*3*2]);
	error = glGetError();
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, 0);

	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, pressure->getTextureID());
	glTexSubImage2D(GL_TEXTURE_RECTANGLE_ARB, 0, x, y, res_m, res_n, GL_RED, GL_HALF_FLOAT_ARB, &data[res_m*res_n*4*2]);
	error = glGetError();
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, 0);

	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, pressure->getRenderTextureID());
	glTexSubImage2D(GL_TEXTURE_RECTANGLE_ARB, 0, x, y, res_m, res_n, GL_RED, GL_HALF_FLOAT_ARB, &data[res_m*res_n*4*2]);
	error = glGetError();
	glBindTexture(GL_TEXTURE_RECTANGLE_ARB, 0);
	glPixelTransferf(GL_ALPHA_SCALE, 1);
}


void FluidGridSimulator::genSamplingList()
{
	switch (pointdisalg)
	{
	case 0:
		{
			//Wang tiles
			glDeleteLists(samplingListId, 1);

			initTime();
			wang.toneScale = m_sampleRatio*20000;
			wang.paintPoints();

			samplingListId = glGenLists(1);
			glNewList(samplingListId, GL_COMPILE);
			glBegin(GL_POINTS);
			int dimMax = vpW; if (dimMax < vpH) dimMax = vpH;
			for (int i = 0; i < wang.numPoints; i++)
			{
				int x = wang.points[i].x * dimMax;
				int y = wang.points[i].y * dimMax;
				if (x < vpW && y < vpH)
					glVertex3f(x + 0.f, y + 0.f, wang.points[i].z);			
			}
			glEnd();
			glEndList();
			printf("%d points (Wang tiles) was generated successfully\n", wang.numPoints);
			break;

		}
	case 1:
		{
/*			//Parallel Poison
			glDeleteLists(samplingListId, 1);
			samplingListId = ps.process(vpW, vpH, 0.001f);
			break;
*/
		}
	case 2:
		{
/*			//Penrose tiling
			glDeleteLists(samplingListId, 1);
			samplingListId = glGenLists(1);
			glNewList(samplingListId, GL_COMPILE);
			glBegin(GL_POINTS);
			PenroseTiling penrose;
			penrose.Tiling(vpW, vpH, m_sampleRatio);
			glEnd();
			glEndList();
			break;
*/

			//halton sequence
			glDeleteLists(samplingListId, 1);
			samplingListId = glGenLists(1);
			glNewList(samplingListId, GL_COMPILE);
			glBegin(GL_POINTS);
			int dimMax = vpW; if (dimMax < vpH) dimMax = vpH; float samplingHaltonN = dimMax * dimMax / (m_sampleRatio);
			densityHaltonN = samplingHaltonN / (dimMax*dimMax);
			int count = 0;
			for (int i = 1; i <= samplingHaltonN; i++)
			{
				int x = haltonsequence(i, 2) * dimMax;
				int y = haltonsequence(i, 3) * dimMax;
				if (x < vpW && y < vpH)
				{
					//float randval = 1 + floor(7.f*(float)rand()/RAND_MAX);
					glVertex3f(x + 0.f, y + 0.f, i * 1.f / samplingHaltonN);			
					count++;
				}
			}
			glEnd();
			glEndList();
			printf("%d points Halton sequence was generated successfully\n", count);
			break;
		}

	}		
	
}